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HE COMBUSTION of gases, in “ite of the apparent 

simplicity, is so complicated a process that many phe- 
nomena that accompany it have not been sufficiently in- 
vestigated. One such phenomenon is the formation of the 
shock discontinuity ahead of a flame front. 

In addition to being of theoretical interest, a study of this 
process also has great practical significance, since the presence 
of a shock discontinuity ahead of a flame front is the condition 
for the transition of a slow combustion wave to a detonation. 
The rate of acceleration of the flame determines where the 
shock wave (and consequently the transition into detona- 
tion) occurs relative to the ignition point in the combustion 
chamber. In the first approximation the problem of the 
motion of a gas ahead of the flame front can be reduced to the 
problem of motion of gas in front of a piston, accelerated in 
accordance with the same law that governs the actual propa- 
gation of the flame. The problem of the combustion of a gas 
contained in an infinite cylinder, bounded on one end by an 
accelerated piston, was first solved by Hugoniot at the end of 
the 19th century. 

In solving this problem, Hugoniot started with the equation 
of motion of an ideal gas, obtained from conservation laws in 
the following form i i 


[1] 


[2] 


Here u is the displacement of a point that has a coordinate x 
at the instant t, y = c,/c,, and a is the velocity of sound in the 
unburned mixture. Disturbances travel in front of the ac- 
celerated piston in the tube, and a discontinuity, namely a 
shock wave, is produced where these disturbances meet each 
other. According to Hugoniot, two elementary disturbances 
produced by the piston at instants that are close to each other 
will meet at a distance 


2a%{1 + [(y — 


X= 3] 

(YF 


t= time of enattatien: of a given disturbance in the gas 
v = velocity of motion of the piston at the instant t 


Hugoniot has shown that all the elementary disturbances 
meet at one point only if the piston moves in a special way, 
— if the velocity of the piston changes with time as 


at\ 2a 
[4] 


In the general case, the disturbances produced by the moving 


Translated from ‘Physical Gasdynamics,’’ USSR Acad. Sci., 
pp. 151-162. 
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piston wil be concentrated not in one section of the tube, but 


over a certain length of the tube. 

It must be noted that by reducing the problem of the mo- 
tion of a gas ahead of the flame front to a problem of the 
motion of a gas ahead of a piston, the flame is assumed to 
have properties of a solid wall, which compresses the gas in 
front of it. 

The time interval within which the elementary disturbance 
sent by the flame evolves into a shock wave is determined by 
the character of the variation of the flame velocity during the 
combustion process. The almost complete lack of experi- 
mental data on the formation of a shock discontinuity ahead 
of a flame front, makes it difficult to check the validity of re- 
ducing the complex problem of motion of gas ahead of a 
flame front to the simpler and already solved problem of the 
motion of gas in front of a piston. 

It should be noted that in most studies of unsteady combus- 
tion or of the transition from slow combustion into detona- 
tion, the process is usually recorded by direct photography of 
the intrinsic radiation of the burning mixture. 

The formation of a shock wave ahead of the flame front is 
not accompanied by intrinsic radiation, and therefore in order 
to follow all the details of this process, it must be examined 
optically, and the visible processes that occur in the unburned 
mixture become decisive. It is best to investigate a hydrogen- 
oxygen mixture, since the combustion process in such mix- 
tures is so rapid, that the shock discontinuity is formed a 
comparatively short distance from the ignition point in the 
combustion chamber. 

To be able to trace all the phases of the process of shock 
wave production, it is necessary to use high speed pho- 
tography. The best known version of high speed photography 
is “spark photography,” in which the subject is illuminated 
by a series of electric sparks produced by discharging a 
capacitor through the illuminating spark gap. 

A distinguishing feature of spark apparatus is the simplicity 
of the spark circuit, the ease of synchronization of the light 
flashes with the investigated process, and the possibility of 
varying the frequency of the photographs over a wide range. 
A spark discharge, lasting 10~ sec, insures a sufficiently short 
exposure time, making it possible to obtain sharp photographs 
even with very fast processes. Inasmuch as the speed of 
propagation of the elementary disturbances that evolve into 
a shock wave does not exceed 1000 m/sec, we have deemed it 
advisable to use, in our particular investigation, a spark de- 
vice capable of taking a total of about 100 frames ata rate of 
50,000 to 70,000 frames per sec. 


1 Experimental Setup 


The experimental setup for studying the process of shock 
discontinuity formation ahead of the flame front consisted 
of an optical section to visualize the disturbances that occur 
in the gas during the combustion of the mixture, an electrical 
section to produce a series of sufficiently intense short dura- 
tion electric flashes, and a photographic section to record 
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the process. The IAB-451 mirror-meniscus apparatus was 
used to make the disturbances in the gas visible. 

The light source used to take pictures by the scanning 
method was a PZh-26 lamp. During high speed spark pho- 
tography, the lamp was replaced by a discharge gap. In the 
latter case the source of light was a series of illuminating 
sparks. The number of sparks in the series and their repeti- 
tion frequency were determined by the parameters of the 
electric circuit. 

The diagram of the electric section of the spark apparatus 
is shown in Fig. 1. The power supply, consisting of a high 
voltage transformer 7’r, which steps up the voltage to 27 kv, 
and of a rectifier Kp, is disconnected prior to taking the pic- 
ture. Thus the source of electric energy is the 1-uf capacitor 
C,. When switch P is closed, capacitor C; is discharged and 
charges capacitor C2, which has a rating of 0.125 uf. When- 
ever the voltage across capacitor C; reaches a value sufficient 
to break down the spark gap, a spark is produced in the gap F. 
Using an inductance L = 2h, a total of 150 to 200 spark flashes 
were produced at a repetition rate of approximately 50,000 
per sec. Kozachenko, using an installation of this type, was 
able to vary the taking rate from 5000 to 100,000 frames per 
sec by varying the inductance from 10 to 0.5h. This spark 
cireuit is convenient in operation. True, a certain difficulty 
is caused by the need of preparing a set of inductors designed 
for rather high voltage. This is why the inductor was later 
replaced by a resistor. If the resistor is chosen sufficiently 
large and C; sufficiently small, the capacitor C; will discharge 
into the illuminating spark gap much more rapidly than it is 
charged through the resistance. This guarantees good quench- 
ing of the spark. By varying the resistance it is possible to 
adjust the picture-taking frequency over a wide range. With 
27 kv across capacitor Ci, and with C,; = 2.7 uf and C, = 
0.125 uf, the breakdown voltage of the gap approximately 15 
kv, and the resistance 1200 ohms, the repetition rate of the 
flashes reaches 50,000 per sec, the total number of flashes being 
approximately 200. The explosion was synchronized with the 
series of illuminating flashes by means of a special pendulum 
switch. 

To obtain bright illuminating flashes, a special discharge 
gap was designed, a so-called “‘linear source of light,’’ shown 
in Fig. 2. The electrodes of the gap are a brass plate 6, 
clamped on both ends by fiber plates 9, and the second elec- 
trode is a tungsten rod 10, which can be moved by means of a 
micrometer screw 11. Light strikes the slit of the optical 
setup through a Plexiglas covered window cut through the 
forward portion of the discharge assembly. Valve 7 allows 
one to pump the air out of the discharge assembly and to fill 
it with hydrogen. In spite of the very low light yield of a 
hydrogen discharge, this gas was chosen, since it permits tak- 
ing pictures at a high frequency. The explosion process was 
recorded on a film having a speed of 250 GOST (= Russian 
standard) units, placed on the outer surface of a uniformly 
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Fig. 2 Discharge assembly, ‘“‘linear source of light?” 


7r 
Fig. 1 Diagram of electrical section of the spark assembly 


rotating drum driven by a synchronous motor. To determine 
the speed of the moving film, an intermittent source of light, 
which produced a series of lines on the film, was photographed 
simultaneously outside the field of the image. The speed of the 
film could be determined from the repetition frequency of the 
light flashes by measuring the distance between the time 
markers. In our case the film moved at a speed of 45 m/sec. 
The pulsating source of light for the time markers was a 
neon lamp MN-7, fed from a 50-500 cycle frequency con- 
verter. To prevent the time markers from becoming super- 
imposed on each other, a shutter was placed in front of the 
neon lamp, which in turn was placed in a special tube; the 
shutter was closed after each revolution. 

The object lens used was an f/1.4, 13.5-cm focal length 
“Sonnar” telephoto lens. It was already noted in the fore- 
going that the process was photographed by two methods: 
A scanning method, and spark motion-picture photography. 
In the former case a continuous source of light was used. To 
prevent superimposed images, a focal-plane type shutter 
was placed in the path of the light and was closed after each 
revolution. The unique construction of the shutter made it 
possible to synchronize the operation of the shutter with the 
process. In the case of spark photography there was no need 
for a shutter. To prevent superimposed frames, the size of 
capacitor C, was chosen such that the entire series of il- 
luminating flashes took place within a single revolution of the 


photorecording drum. 


2 Explosion Chamber and Sequence of A 
Experiments 


To explain the effect of the shape of the chamber cross sec- 
tion on the process of formation of a shock wave in front of the 
flame front, the behavior of the process in round and rec- 
tangular chambers was compared. The area of the transverse 
cross section of both chambers was kept constant. Cham- 
bers in which the combustible mixture was burned consisted 
of several sections joined with flanges. By using a different 
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number of sections, it was possible to vary the length of the 
explosion chamber from 280 to 608 mm. The chamber was 
closed on both ends with cover plates. A spark plug to ig- 
nite the mixture was mounted in one of the plates, and a valve, 
through which the air was evacuated and the chamber filled 


20 mm 120 mm 


Fig. 3 Instantaneous photographs of the combustion of hydro- 
gen-oxygen mixtures in the first section of a rectangular chamber 


Fig. 4 Formation of shock discontinuity ahead of the flame 
front 


with the explosive mixture, was mounted in the other plate. 
The process was photographed through side openings cut in 
one of the sections of the chamber. 

The openings were covered with plane-parallel optical-glass 
plates, secured to the body of the chamber with the aid of 
slotted covers. The cross pieces of the covers served as scale 
markers. The combustible mixture was compounded by 
volume and stored in a gas holder without being dried before 
filling the chamber. We investigated the combustion of 
hydrogen-oxygen mixtures containing 25, 33.3, 50, 66.6 and 
70 per cent hydrogen. 

After the chamber was filled with the explosive mixture, the 
electric section of the apparatus was turned on. The ca- 
pacitor bank was charged to 27 kv, after which the motor of 
the drum camera was turned on and the latch, which main- 
tained the pendulum of the synchronizing device in a horizon- 
tal position, was actuated. The falling pendulum switch 
closes the high voltage circuit of the spark assembly. Simul- 
taneously, the tail piece of the pendulum opens the circuit 
of the primary winding of the transformer, producing a spark 
that ignites the mixture. Approximately 200 different phases 
of the investigated process were recorded on the uniformly 
moving film of the photorecorder. 

The first series of experiments was carried out in chambers 
long enough to permit a shock wave to be formed before the 
disturbances, produced by the front of the flame in the un- 
burned gas, reached the end of the chamber. 

The second series of experiments was carried out in a short 
chamber. In this case a shock wave is formed as a result of 
cumulation of disturbances reflected from the end of the 
chamber. 


3 Description of the Experimental Material and 
Analysis of the Results 


a Formation of Shock Wave Ahead of the Flame Front 


Fig. 3 shows a series of instantaneous photographs of the 
combustion of hydrogen-oxygen mixtures in the first section 
of a rectangular chamber measuring 19 X 37 mm. 

The photographs show that the flame propagates from left 
to right as the film moves in a vertical direction. The two 
black strips parallel to the time axis are scale markers, spaced 
38mm apart. The photographs were taken against a so-called 
“gray” background. 

As can be seen from the photographs, the initially hemi- 
spherical flame front is gradually drawn out and becomes 
meniscus shaped. As the flame propagates, the flame front 
recorded on the photograph becomes somewhat thicker, show- 
ing development of a flame front surface with the meniscus 
shape retained. 

The flame front is preceded by traveling disturbances, which 
are quite clearly seen on photographs of Fig. 4, which charac- 
terize the development of the process of combustion of a hy- 
drogen-oxygen mixture containing 70 per cent hydrogen in 
the second section of a chamber 480 mm long. The distance 
from the point of ignition to the start of the section viewed 
was 160 mm. The developing shock wave appears as a 
gradually narrowing bright zone, which changes finally into 
a thin line, evidencing the formation of a shock discon- 
tinuity. The shock discontinuity is formed approximately 
25 cm from the spark plug that ignites the mixture. In a 
mixture containing 50 per cent hydrogen, the shock discon- 
tinuity is formed at a distance of 29 to 30 cm, and in a mix- 
ture with a still greater oxygen content, the elementary dis- 
turbances sent by the flame evolve into a shock wave 57 to 
58 cm from the spark plug. If the flame front traveling 
through the tube can be likened to a piston that compresses 
the gas in front of it, then the already available Hugoniot 
solution can be used to determine the distance from the 
ignition to the place of formation of the shock discontinuity. 
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To estimate the value of X in Equation [3] for the case of 
propagation of the flame in a tube, it is necessary to know the 
speed and acceleration of the flame, the + of the mixture, and 
the velocity of sound propagation in the unburned mixture. 
The velocity of sound is calculated from the formula - ri 


a = V981y7RenT 


where 
T = absolute temperature of the mixture 
Rem = Ro,*Ru,/(aRo, + BRx,) 
Ro, = 26.5 m/deg 

Rs, = 420 m/deg 

a fraction of the hydrogen in the mixture ; 
B fraction of oxygen in the mixture a 


For the hydrogen-oxygen mixtures of interest to us, the 
velocity of sound at 20 C has the following values: 


in mixture 4, m/sec 

70 (567 

535 

B94 

The velocity and acceleration of the flame at any instant 
of time can be readily found experimentally. During the first 
instants after the ignition of the mixture, the speed of the 
flame first increases and then remains constant, after which it 
gradually decreases and drops sharply at the instant when 
the flame makes contact with the walls of the chamber. This 
is followed by a new increase in the speed of flame propagation. 
In a rectangular chamber, by virtue of the fact that the flame 
front first touches the upper and lower walls and then the 
side walls of the chamber, it is impossible to establish a reduc- 
tion in the speed of flame propagation due to the increase in 
heat transfer during contact between the flame and the 
walls. In later instants of time the speed of flame propagation 
increases as the flame advances in the tube, in both round and 
rectangular chambers. If one disregards the stage during 
which the flame propagates from ignition to contact with the 
wall, the path covered by the flame in a time ¢ can be repre- 
sented by the empirical formula 


= path, m 
t = time, sec 
a,b,c = quantities that depend on the concentration of the 


mixture and structural features of the chamber 


Differentiating Equation [5] with respect to ¢, we can 
readily find the velocity and acceleration of the flame at any 
instant of time. Figs. 5a and 5b illustrate the acceleration as 
a function of time for hydrogen-oxygen flames in round and 
rectangular chambers. 

It is easy to note that the acceleration of the flame depends 
on the structural features of the chamber, and depends for 
each given chamber on the concentration of the mixture. 
Whereas in a rectangular chamber the combustion process 
proceeds in such a way that the acceleration of the flame di- 
minishes with time for all the hydrogen-oxygen mixtures con- 
sidered, in a round chamber of equal transverse area as the 
rectangular chamber, the dependence of the acceleration of 
the flame on the concentration varies quite pronouncedly. 
In fact, in a mixture containing 70 per cent hydrogen, the 
acceleration of the flame increases with time; in a mixture 
containing 50 per cent hydrogen, the flame propagates with 
equal acceleration, and in a mixture with a still greater con- 
tent of oxygen, the acceleration diminishes gradually as the 
flame moves in the chamber. Since the variation of the ac- 
celeration of the flame should be greatly affected by the dis- 
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Fig. 5 Change in acceleration of flame with time = 


tance from the ignition point to the place where the shock 
discontinuity is generated, it can be deduced without calcula- 
tion that the conditions in a round chamber are obviously 
more favorable for the occurrence of a shock discontinuity 
than in a rectangular chamber. By way of example let us 
calculate at what distance from the ignition point will there 
be encountered the disturbances propagated by the flame 
front during the instant ¢ = 10~‘ sec for the case of a com- 
bustion of a mixture containing 70 per cent hydrogen in a 
round chamber. 

Inserting the corresponding values of v and ¥ in the formula 
for X, and taking into account the fact that the velocity of 
sound in the unburned mixture is 567 m/sec, we obtain X = 
22.5 cm. Considering that at the instant ¢ = 10~‘ sec the 
front of the flame is at a distance 2.4 cm from the edge of the 
chamber, the distance from the point of ignition to the point 
of “generation” of the shock discontinuity is 24.9 cm. This 
value is in good agreement with experiment, thus confirming 
the correctness of the idealization we have assumed for the 
process. The solution of this problem becomes particularly 
clear if the graphic method of characteristics is employed. 


b_ Use of the Graphic Method of Characteristics in the 
Solution of the Problem of Unsteady One-Dimensional 
Flow of Gas Ahead of the Flame Front 


Without dwelling on the derivation of the equation of 
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characteristics, which is covered in detail in (2 and 3), let us 
proceed to solve the problem of interest to us. 

A flame front propagates in a hydrogen-oxygen mixture 
contained in a tube at constant pressure and temperature. 
The flame propagation is described by the curve OA in the 
(x, #) plane. It is required to find the state of the gas ahead 
of the flame front. Let us make the corresponding graphic 
construction. The sonic disturbances at the instants of time 
1, 2, 3, 4, 5, 6, 7 and 8 are determined in the (v, a) plane by 
the characteristic represented in Fig. 6c by a straight line 
making an angle 45 deg with the v and a axes. This angle is 
determined by the choice of the scale of the ordinate axis. 

The sonic disturbances traveling from points 1, 2, 3,...,8 
(Fig. 6a) are determined as the ordinates of the points, the 
abscissas of which are the corresponding values of the velocity 
of the flame front. 

Let us call attention to one peculiarity in our problem. 
Usually, in the solution of the problem of motion of a gas in 
front of a piston moving in a tube, it is the first disturbance, 
produced by the piston at the initial instant of time and 
propagating with the velocity of sound in the unburned gas, 
which plays the major role. In the ignition of a mixture by 
even a very weak spark, a disturbance is also produced in the 
gas and propagates at the velocity of sound, but this is where 
the analogy ends. Experience shows that the first disturb- 
ance, produced in the gas by an ignition spark, usually does 
not participate in the formation of a sufficiently intense shock 
wave. Therefore, in solving the problem by the method of 
characteristics, this disturbance should not be taken into 
account. The results of the calculations of the sonic dis- 
turbances, produced by the front of the flame in later instants 
of time, are listed in Table 1. 

As can be seen from the diagram, the characteristics that 
emerge at the instants of time 1 and 2 are parallel, since the 
flame velocity in this section remains constant. The flame 
velocity decreases then, causing a diverging bundle of charac- 
teristics. Starting with the instant t = 10~‘ sec, the flame 
velocity increases. The characteristics that emerge at the in- 
stants 4, 5, 6, 7 and 8 intersect at a distance ~24.5cm. The 
calculation is in good agreement with experiment. Analogous 
calculations were made for combustion of a mixture containing 
50 and 33.3 per cent hydrogen. 

Calculations carried out by the same method show that 
when a flame propagates in a rectangular chamber, the shock 
discontinuity is formed at a considerably greater distance 
from the point of ignition than if the same mixture is ignited 
in a round chamber. 

The results of the construction made by us are presented in 
Table 2. 

It is easy to note that the distance from the point of ig- 
nition to the point of formation of the shock discontinuity de- 
pends substantially on the structural features of the chamber. 

For equal transverse areas, a shock discontinuity is pro- 
duced in the round chamber at half the distance from the 
point of ignition than when the same mixture is ignited in a 
rectangular chamber. 

To relate the distance from the point of ignition to the point 
of production of a shock wave and the time of its occurrence 
with the fuel characteristics of the mixture, we make use of 
the formulas obtained by Hugoniot to calculate the values of 
X and 7. According to Hugoniot 


a 


It has been shown previously that the acceleration of the 
flame, which enters into the denominators of both expres- 
sions, is a function of the construction of the chamber, and is 
determined for a given chamber by the concentration of the 
combustible mixture. 


1 Numbers in parentheses indicate References at end of paper. 


y¥-1 v+a, 

Point v, m/sec m/sec a, m/sec m/sec 
1 3135 627 927 
25 3135 627 927 
3 3035 607 807 
4 2976 595 736 
5 3022 604 791 
6 3070 614 849 
7 3117 623 905 
8 3169 634 968 


Calculation shows that during that time when the front 
of the flame generates disturbances in the unburned mixture 
that cumulate in a shock wave, the acceleration of the flame, 
in the range of concentrations of interest to us, is propor- 
tional to the square of the percentage of hydrogen in the 
mixture and, consequently 


Fig. 7 shows the dependence of X on a*/m*. The abscissas 
represent the square of the ratio of the velocity of sound to 
the percentage of hydrogen, while the ordinates represent 
the distance from the point of ignition to the point of genera- 
tion of the shock discontinuity. The open circles denote the 
values of X that correspond to the process in a round cham- 
ber, while the filled ones are for a rectangular chamber. 

So far we have not considered the process of formation of 
the shock wave in chambers that are not sufficiently long. 
In short chambers the elementary disturbances, generated by 
an accelerated flame, do not have a chance to cumulate into a 
shock wave before they intersect the end of the chamber. 
The formation of the shock discontinuity in a short chamber 
is therefore complicated by the participation of disturbances 
reflected from the end of the chamber. 


ec Formation of Shock Wave in a Short Chamber 


The shock wave is formed at a definite length of the cham- 
ber. If the disturbances produced by the flame in the un- 
burned gas do not have time to cumulate into a shock wave 
prior to encountering the end of the chamber, the shock dis- 
continuity may occur as a result of addition of the dis- 
turbances that are reflected from the end of the chamber or 
from the front of the flame. Elementary disturbances, as a 
rule, do not pass through the front of the flame, and are re- 


Rectangular 
mere, Round chamber————. 
aor Distance from Distance from Distance from 
Bae, ignition point ignition point ignition point 


to point of 
formation of 


to point of 
formation of 


to point of 
formation of 


Per cent shock wave, shock wave, shock wave, 
H, in cm cm cm 
mixture (calculated) (experimental) (calculated) 
70 24.5 25-26 
50 28.5 29-30 69 
33.3 56 Po =58 121 
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Fig. 6 Solution of the problem of unsteady one-dimensional 
flow of gas ahead of the flame front by the method of charac- 
teristics 


LT 


Fig. 7 Dependence of the value of X on the square of the ratio 
of the velocity of sound to the percentage of hydrogen 


flected from it as from a solid wall. Fig. 8 shows the process 
of formation of a shock discontinuity in the combustion of a 
mixture containing 33 per cent H; in a chamber 280 mm long. 
In the left half of the diagram is shown the trace of motion of 
the flame front in the first section of the chamber. The de- 
veloping shock wave is clearly seen ahead of the flame front 
that appears in the field of view on the twelfth frame. On the 
first frames it has the form of a rather broad light zone. After 
reflection from the end of the chamber the disturbance propa- 
gates in the opposite direction and is recorded on the photo- 
graph in the form of a gradually narrowing dark zone. In 
addition to the formed shock wave, the photograph shows 
clearly the disturbances produced by the flame front in the 
unburned gas. Propagating between the flame front and the 
end of the chamber, they form a dense grid of disturbances 
that do not interact with each other. 

Since the flame front behaves like a solid wall with respect to 
the elementary disturbances, it is possible, when calculating 
the state of the gas between the front of a flame and the wall, 
JaNuaRY 1960 BA, 


Fig. 8 Formation of shock discontinuity in a short chamber 


to solve the problem by the graphic method of characteris- 
tics, by likening the flame to a piston. Fig. 9 shows the 
corresponding constructions for the case represented in the 
photograph. The law of propagation of the flame is specified 
by curve OA. The position of the wall is noted by a vertical 
line with coordinate X = 25mm. Comparison of Figs. 8 and 
9a shows that the calculation of the place of formation of 
the shock wave in a short chamber is in good agreement with 
experiment. 
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Conclusions 


1 A flame moving with accelerated motion along a tube 
produces disturbances in the unburned gas. A series of such 
disturbances, cumulating in one section of the tube or over a 
certain length of it, produces a shock wave, clearly seen on the 
photographs as a gradually narrowing bright zone, which 
finally becomes a thin line. 

2 The formation of the shock discontinuity occurs over a 
definite length of the chamber. If the disturbances produced 
by the flame in the unburned gas do not have a chance to 
cumulate into a shock wave prior to encountering the end of 
the chamber, the shock discontinuity may occur as a result 
of the addition of disturbances reflected from the end of the 
chamber or from the front of the flame. Elementary dis- 
turbances, as a rule, do not pass through the flame front, but 
are reflected from it as from a solid wall. 

3 In the calculation of the process of formation of the 
shock wave ahead of the flame front it is possible, by reducing 
the problem of the motion of gas in front of the flame to the 
problem of motion of gas in front of a piston accelerating in 
accordance with the same law as the flame actually ac- 
celerates, to make use of the Hugoniot formulas to deter- 
mine the distance from the point of ignition to the point of 
formation of the shock discontinuity, or else to solve the 
problem graphically by the method of characteristics. Our 
own calculations were in good agreement with experiment, 
this being evidence of the correctness of the idealization as- 
sumed for the process. 

4 The dependence of the distance X from the point of ig- 
nition to the point of formation of the shock discontinuity 
on the concentration of the mixture, can be described, in the 
region of concentrations investigated by us, by the relation 
X ~ a*/m?, where a is the velocity of sound in the mixture 
and m is the percentage of hydrogen in the mixture. 

5 The distance from the point of ignition to the point of 
formation of the shock discontinuity depends on the shape of 
the cross section of the chamber. In a round chamber the 
acceleration of the flame front is much greater than in a 
rectangular chamber (for equal transverse areas). By virtue 
of this, the value of X in the former case is almost 2.5 times 
smaller than in the latter. 
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is accompanied by a wave process. The acoustic waves 
are the cause of one of the unique propagations of the flame, 
namely, vibrational combustion. The ripples of the flame are 
connected with phenomena of hydrodynamic character that 
occur during combustion, which are much more intense than 
pure acoustic phenomena. The first attempt to attribute 
the occurrence of flame ripples to purely mechanic interaction 
of the disturbances that arise in the gas with the front of the 
flame was made by Dixon, who attributed the flame ripples to 
the interaction between the reflected sound wave and the 
front of the flame. By using Dixon’s method to plot the 
path of the sound wave, Sokolik found no agreement between 
the reflected wave and the place where the flame is rippled, 
which is quite understandable, since the sound wave is a wave 
of infinitesimally small amplitude and cannot cause finite 
ripples of the flame. Flame ripples are due to the action of 
shock waves on the front of the flame. In spite of the fact 
that the combustion with ripples has long interested the 
researchers, the problem of the investigation of the interaction 
of the flame front with the shock wave is still not being given 
enough attention at the present time. Most investigations 
devoted to this problem pertain to a single particular case, 
which indeed is of great practical and theoretical interest, 
namely the case where the interaction of the shock wave with 
the flame of the front results in a detonation wave. Com- 
pletely opposed to this is the problem of interaction between 
a shock wave and the flame front without detonation propaga- 
tion of the flame. Without dwelling in detail on the mecha- 
nism of formation of shock waves ahead of the flame front, 
considered elsewhere,! we note that since the shock dis- 
continuity occurs at a definite position in the chamber, the 
front of the flame can interact both with the shock wave that 
is being formed and with shock waves that have already been 
formed. The interaction process will have, in either case, its 
own characteristic features. To follow the details of this 
process, it is necessary to visualize the processes occurring in 
a fresh mixture and not accompanied by intrinsic glow. 
During the process of interaction, changes occur both in the 
flame front and in the shock wave, and to record the process 
it is necessary to use high speed motion picture photography, 
which makes it possible to examine all its phases. 

A detailed description of the optical portion of the equip- 
ment with which one can make visible even weak disturbances 
produced by the flame in a fresh mixture has been given.} 
The same article describes the high speed spark apparatus, 
which yields up to 200 different phases of the investigated 
process at a rate of 50,000 frames per sec. 

We used the apparatus to investigate the interaction be- 
tween the front of a flame propagating in closed chambers and 
a shock wave produced ahead of the flame front as a result 
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of said of iii produced by the flame in the 
fresh gas. Use was made of hydrogen-oxygen mixtures, 
containing 25, 33.3, 50, 66.6 and 70 per cent hydrogen. The 
combustion chamber length varied from 280 to 480 mm. The 
pictures were taken through windows cut in the side walls of 
rectangular chambers 19 X 37 mm. The windows of the 
chambers were covered with optical glass. The metal strips 
used to clamp the glass to the body of the chamber were used 
as scale markers, which were fixed on the photographs in the 
form of dark strips, parallel to the time axis. The process was 
photographed by two methods: The method of scanning, and 
the method of high speed motion picture photography. 

After these preliminary remarks, let us proceed to examine 
the problem of the interaction between the flame front and 
the shock discontinuity. 


1 Interaction Between the Flame Front and 
the Shock Discontinuity 


Fig. 1 shows photographs of the process of combustion 
of hydrogen-oxygen mixtures, containing 50 per cent hydrogen 
in a chamber 280 mm long, obtained by the scanning method. 
Recorded to the left are the development of the process in the 
first section of a two-section chamber, and to the right are the 
processes near the second end surface. The flame accelerated 
in the tube excites in the fresh mixture disturbances which 
cumulate, after reflection from the end of the chamber, into a 
shock wave which passes freely through the front of the flame, 
almost stopping the latter. The flame, so to speak, “splashes 
away.” The second shock wave, formed as a result of cumu- 
lation of disturbances propagated by the flame in later in- 
stants of time, will be reflected by the end of the chamber and 
will encounter the already “splashed”’ flame, through which 
this shock wave also passes freely, pushing its front back some- 
what. Photographs obtained by the scanning method dis- 
close nothing concerning the changes that occur in the flame 
front during the process of this interaction. Fig. 2 shows 
a series of photographs of the same process, obtained by spark 
photography at 50,000 frames per sec. The disturbances that 
evolve into a shock wave are seen distinctly ahead of the 
flame front, which has a clearly pronounced meniscus shaped 
form. 

The encounter between the shock wave and the front of the 
flame is accompanied not only by a reduction in the speed of 
propagation of the flame, but also by a deformation of the 
latter. A small indention appears in the heretofore meniscus 
shaped front of the flame, which now assumes a characteristic 
tulip-like shape. After the passage of the shock wave through 
the front of the flame, a long comet-like loop stretches out 
behind the flame front, in which the mixture is afterburned. 
The unique shape of the “loop”’ is evidence of a definite dis- 
tribution of the velocity field of the stream, of the flow, ac- 
companying the shock wave. The front of the flame acquires 
a cellular structure. The total surface of the front increases 
and this, in turn, increases the rate of propagation 
of the flame. Gradually the flame front again acquires a 
meniscus shape. The second shock wave, encountering the 
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Fig. 1 Scanning of the combustion of a hydrogen-oxygen mixture 


front of the flame and accompanying the comet-like “loop,” 
is in itself deformed somewhat after passing through the region 
in which the mixture is afterburned. Frames 35 to 37 of 
photograph 2b shows clearly the curving of the front of the 
shock wave, following the comet-like “loop.’”’ Upon leaving 
this zone, the front of the shock wave again becomes flat. 
Since the shock wave traveling ahead of the front encounters 
the flame several times before the latter reaches the end of the 
chamber, the process of variation of the shape of the front of 
the flame during combustion with ripples is periodic in charac- 
ter. 

We considered above the interaction between the flame 
front and the shock wave, producing the so-called combustion 
with ripples. At a sufficient length of the combustion cham- 
ber, the flame front may catch up with the shock wave before 
the latter reaches the end of the chamber. Under certain 
conditions imposed on the intensity of the shock wave and on 
the speed of the reaction inside the flame front, such an inter- 
action between the wave and the flame front results in the 
formation of a detonation wave. The transition into a detona- 
tion wave is observed also during combustion of the mixture 
in a chamber, the length of which is less than the length of 
the predetonation path. For this transition it is essential 
that there exists an intense shock wave and that the mean 
velocity of the reaction be capable of maintaining the propa- 
gation of the detonation front. 

The existence of an intense shock wave is a necessary but 
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Fig. 2 Series of spark photographs of the same process 


not sufficient condition for the transition from slow combus- 
tion into detonation, as evidenced by the series of photo- 
graphs in Figs. 3a and b, which record the combustion of a 
mixture containing 50 per cent hydrogen near the end surface 
of a chamber 480 mm long. In spite of the fact that in both 
cases a shock wave of the same intensity travels in opposition 
to the front of the flame, in the former case the process is such 
that the encounter between the shock wave and the flame 
front does not lead to the formation of a detonation wave. 
In the second case the slow combustion does change into 
detonation. Unfortunately, the photographs did not register 
the interaction between the wave and the flame front, since 
the latter took place in a section of the chamber obscured 
from view. 


2 Interaction Between the Forming Shock hdd 
Wave and the Flame Front 
If the shock wave passes freely through the flame front, the 


latter behaves like a deformable wall with respect to the form- 
ing shock wave. The process of the interaction between the 


forming shock wave and the flame front can be followed on the 
series of photographs (Fig. 4) which show the combustion of a 
mixture containing 33.3 per cent hydrogen in a chamber 280 
mm long. The forming shock wave, traveling opposite to the 
flame front, can be seen in the form of a gradually narrowing 
dark zone. 


After encountering the flame front, it deforms 


ARS JOURNAL SUPPLEMENT 


i 
> 
» 
| 


120 


120 


> 


Fig. 3 Interaction between flame front and shock wave — 


the latter and cannot pass through it. The photographs show 
clearly how a meniscus shaped flame front becomes tulip 
shaped under the influence of a shock wave. The deforma- 
tion of the flame front when a wave is reflected from it is 
probably the cause of the change of the angle of inclination 
of the front of the forming shock wave to the chamber axis 
after the reflection. 


The interaction between the forming shock wave and the 
front of the flame can lead to the formation of a new ignition 
center ahead of the flame front. This process is quite notice- 
able on the series of photographs (Fig. 5a) showing the com- 
bustion of a mixture containing 25 per cent hydrogen in the 
second section of a 280-mm long chamber. The new ignition 
center, the appearance of which is quite noticeable starting 
with the fourth frame, has the form of a gradually broadening 
dark region. This is followed by a bright region—the com- 
pression region, contained between two combustion centers. 
The occurrance of the new center is due apparently not to the 
ignition of the mixture by compression in the wave, but by the 
injection of the reaction products in the fresh mixture, as 
evidenced by the thin bridge between the new ignition center 
with the flame front. 

As can be seen from the photographs, the production of a 
new ignition center does not change substantially the charac- 
ter of the combustion. The interaction between the second 
forming shock wave and the flame front leads, as can be seen 
from the last frames of the photograph 5, to a new change in 
the form of the front and to a sharp increase of its surface, 
owing to the appearance of the cellular structure. The 
cellular structure is particularly clearly pronounced on photo- 
graphs (Fig. 5b) of the combustion of the same mixture in 
the second section of the three-section chamber, 480 mm long. 

An examination of the above experimental data leads to the 
following conclusions. 

1 The forming shock discontinuity, traveling in opposition 
to the flame, does not as a rule pass through the flame, but is 
reflected from it. The encounter between the forming shock 


Fig. 4 Interaction between flame front and forming shock wave 


Fig. 5 Formation of a new ignition center: a—ahead of the flame front; b—cellular structure of flame front 
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wave and the flame front is accompanied by a deformation of 
the latter. In some cases a new ignition center is formed 
ahead of the flame front, produced by the injection of the 
reaction products into the fresh gas. The shock wave itself 
also experiences a change. 

2 The produced shock wave can readily pass through the 
flame front and throw it back. The throwback of the flame is 
accompanied by a deformation of its front. In a flame front 
which is meniscus shaped before the encounter with the shock 
wave, there is produced a dent which imparts to it a charac- 
teristic tulip shaped form. The structure of the flame front 


PAPER by B Bolz published i in 1955, is 
devoted to the influence of turbulence on the speed of 
flame propagation. It can be shown that the paper contains a 
principal error in the processing of the experimental data, an 
error leading to a considerable underestimate of the turbu- 
lent speed of combustion propagation; furthermore, the ex- 
tent of the underestimate is different in different cases. 
Since reference may be made to this paper when comparing 
various turbulent combustion theories with experiment, we 
deem it advisable to analyze this paper in greater detail. 
Bolz and Burlage, igniting a jet of homogeneous mixture of 
commercial methane with air by means of a single spark, 
produced a small spherical flame, carried by a stream having 
an equalized velocity profile. By taking successive photo- 
graphs of the growing spherical flame, they determined the 
apparent combustion propagation speeds in the absence of 
walls. The speeds of the flame relative to the uncombusted 
gas wu were calculated by dividing the apparent speed ug by 
the increase in volume in combustion at constant pressure C 


where mp is the number of moles, Rp the gas constant, 7) the 
absolute temperature of the uncombusted gas, and n, R and T 
are the same quantities for the adiabatically expanded com- 
bustion products. Under the conditions of the Bolz and 
Burlage experiments, the ratio of the volume of the com- 
busted gases to the initial volume of the combustible mixture 
C ranged from 6.65 to 7.47, depending on the composition of 
the mixture. 

The speeds of laminar flame for various mixture composi- 
tions, determined by Bolz and Burlage, were in good agreement 
with data obtained by others. The values of the speeds of 


propagation of turbulent combustion were not expected by __ 


Translated from Izvestiia Akademii Nauk SSSR, Otdelenie 
Tekhnicheskikh Nauk, Energetika i Avtomatika (Bull. USSR Acad. 
Tech. Sci., Div., Power and Automation), no. 2, 1959 ,pp. 137- 
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of Combustion 


is also changed, for the front becomes covered with cells that 
increase considerably the combustion surface. 

A long comet-like “loop” extends behind the flame front, 
in which the mixture is afterburned. The increase in the sur- 
face of the flame leads to an increase in its propagation speed. 
The new encounter between the shock wave and the flame 
causes a new deformation of the flame front, etc., making it 
possible to speak of a periodic change in the shape of the 
flame front in combustion with ripples. In passing through 
the flame front and the comet-like loop behind it, the front of 
the shock wave becomes curved, and resumes the initial form 
after emerging from the zone of the “loop.” 
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the ‘idles In some cases, en in mixtures with 
very lean or very rich combustible components, the turbulent 
speeds were found to be unusually low, even less than normal 
speeds of flame propagation in the same mixtures. The good 
reproducibility and sufficient compactness of the results made 
it impossible to attribute this unexpected result to random 
measurement errors. 

We shall show that the Bolz and Burlage experiments con- 
tain a systematic error in the processing of the experimental 
data, due to the neglect of the width of the combustion zone. 

Equation [1] is valid only when the width of the combustion 
zone is small compared with the radius of spherical flame. 
The combustion of the gas is still incomplete in the combus- 
tion zone, and if its width is comparable with the dimensions 
of the flame, the expansion of the flame is not determined by 
the quantity C, which ranges from 6.65 to 7.47, as is suggested 
by Bolz and Burlage, who believed that the spherical flame 
contains the products of complete combustion. It is deter- 
mined instead by a smaller quantity; furthermore, the wider 
the combustion zone, and the more uncombusted gas it con- 
tains, the smaller this quantity. 

The width of the combustion zone in the Bolz and Burlage 
experiments on laminar flames was always small compared 
with the flame dimensions. It was on the order of 0.2 mm, 
while the flame measured several dozens of millimeters. 
However, when measuring the turbulent speeds of flame 
propagation, the foregoing conditions were not satisfied—the 
width of the turbulent combustion zone amounting to a notice- 
able fraction of the flame radius. 

For the case when the scale of turbulence 1 exceeds the 
width of the laminar combustion zone, the size of the turbu- 
lent combustion zone \ can be determined semi-empirically 
from the expression 


= BU v/tin)™ [2] 


Here v/up is the ratio of the turbulent pulsational compo- 
nent of flow velocity to the normal flame velocity, and B is a 
constant ranging from 1 to 3.5. The exponent m assumes 
values between 1 and 0.5, depending on the theory of turbu- 
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lent combustion employed; in the experiments described here, 
m is probably close to 0.7 or 0.8. 

Strictly speaking, Equation [2] applies only when v >> wa, 
while in this case v is commensurate with u,. However, for 
lack of anything better, we shall use Equation [2] to estimate 
the extent of the turbulent combustion zone, making use of 
the values of v and uz listed in Table 4 and values of / from 
Table 5 of the cited paper. Setting B equal to the average 
value, 2.25, we obtain for experiments 63 and 66 (for which 
lowest speed of Table 4 is given) the following values for the 
width of the zone of turbulent flame 


16.4 0.76 
' = 2.2. 2. = 5. 
Aes 5.0 mm 
16.4 
= 2. 
= 2.25 X 


These values cannot be neglected compared with the flame 
dimensions, if we recall that the maximum flame radius does 
not exceed 30 mm. 

For an approximate estimate of the influence of the width 
of the zone on the quantity C, we can assume the following 
simplified scheme for the structure of the front: Half the 
combustion zone is filled with combustion products, while the 
other half contains uncombusted gas. Under these assump- 
tions, the coefficient C,, corrected for the zone width, is re- 
lated to the coefficient C used by Bolz and Burlage, as follows 


Cc 
(1 — + — (1 — 


Assuming that the measurements were made at maximum 
flame radius ? = 30 mm (see Fig. 7 of the cited paper), taking 
the values of C from Table 4, and using the values of \ as com- 
puted above, we obtain for experiments 63 and 66 the values 
of 0.45 and 0.4 respectively, for the ratio of the corrected to 
uncorrected coefficient, C,/C. Taking these values into 
account, we find that in experiment 63 the speed of turbulent 
flame propagation was not 24 cm/sec, as indicated in Table 4, 
but 24/0.42 = 56 cm/sec, whereas the normal velocity in the 
same mixture is 21.6 cm/sec. In the case of experiment 66, 
the correction for the width of the zone increases the value of 
turbulent velocity from 18 cm/sec listed in the table to 45 
cm/sec, the normal velocity in the same mixture being 17.4 
cm/sec. Analogous corrections for the width of the zone of 
turbulent combustion should be made for all the experimental 
points. 

It is interesting that the dependence of the flame velocity on 
its dimensions is clearly seen on the Bolz and Burlage curves 
in Figs. 7, 14,15 and 16. As the radius of the flame increases, 
the radius vs. time curves bend as if the speed of the flame 
were to increase with time, while actually one should expect a 
decrease in speed, owing to the damping of the turbulence. 
Bearing in mind that the flame dimensions were too small for 
self-turbulization of the combustion front to arise, the unusual 


[3] 
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= 5.8 mm vozsl¥. 


course of the curves can be taken to be indirect proof of the 
influence of the width of the combustion zone on the expan- 
sion of the combustion products: As the flame dimensions 
increase, the value of C, increases (see Eq. [3]). At a suffi- 
ciently large flame radius, C, goes into C and the apparent 
speed approaches the asymptotic value—the radius vs. time 
curve approaches a straight line. 

Actually, the corrections for the width of the combustion 
zone can hardly be accurate, since the width of the combustion 
zone was not found in the measurements. But even if it were 
known, the correction would be only approximate, since the 
distribution of the combustion products over the turbulent 
flame zone is unknown. The only accurate way of deter- 
mining the turbulent flame speed in free flow, by the pro- 
cedure used by Bolz and Burlage, would be to make the 
measurements at sufficiently large flame sizes, so that the 
width of the combustion zone could be neglected by com- 
parison. 

The remarks made in this communication apply also to 
measurements of turbulent speed in an inverted cone of turbu- 
lent flame. 

In our opinion the most accurate measurements, although 
not free of shortcomings, are measurements in a straight flat 
cone of a turbulent flame. However, in this case one must 
estimate the speed of propagation of turbulent combustion 
from the dimensions of the surface that envelops the forward 
edge of the flame, rather than from the area of the maximum 
glow zone, as was done, for example, by Summerfield, Reiter, 
Kebely and Mascolo (2). To explain the last remark, let us 
perform an imaginary experiment. Imagine a straight flame 
cone. We maintain in it a constant speed of propagation of 
turbulent flame, changing if necessary to pulsational com- 
ponent of the gas flow. To satisfy the condition of invariance 
of the turbulent flame speed propagation, we increase sharply 
the width of the combustion zone, say by increasing the scale 
of turbulence or by reducing the normal flame speed. As the 
width of the zone of turbulent combustion is increased, all the 
characteristic dimensions of the combustion zone increase, the 
upper zone of maximum glow moves away from the mouth 
of the burner and the surface of this zone increases. Using 
the procedure of Summerfield et al. we connect the increase in 
the surface of the zone of maximum glow with the reduction 
in the speed of propagation of turbulent combustion. Actu- 
ally, however, there was no reduction whatever in the speed of 
the flame, and the increase in the zone of maximum glow is 
due to the increase in the width of the combustion zone, 
determined by factors other than the — of propagation of 


turbulent flame. aan 
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Calculation of Orthotropic Conical | = 
Shell for Arbitrary External Load, | = 
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:, _ The conical shell is the component part of most aviation structures. Therefore the knowledge 
Sd reise Ciencias: of an elastically loaded and deformed state of such a shell is very important for the solution of many 


aad niga eels problems. Thus far most papers have been devoted to the strength of conical shells under axially 
Paris =e ees. symmetricalloading. Thecase of arbitrary loading has not been sufficiently considered. The pur- 
=) es Pao Nute pose of this paper is to obtain equations for the design of a conical shell, a type of a rocket nose, 
ws having nonaxially symmetrical loading. 
E SHALL determine the stressed and deformed state of these assumptions become 
an orthotropic conical shell, using the general technical co 
theory of Vlasov (1).! Starting with the physical hypotheses 2. (rN) — kN: + 
of this theory pertaining to orthotropic semi-momentless shells dz 26 a4 
of medium length, we shall consider only the transverse bend- an oM, 
ing moments from among all the internal moments. A— +2 — (rS) +kS +Q,+ ArY = — + 1Q, =0 
The equations of equilibrium of an element of shell under 00 28 
+ 2 + kM, — ArQ, = 0 
where 
2, 0 = cylindrical coordinates of the shell (see Fig. 1) 
r = kz, radius of conical shell 
k = tan ¢, tangent of the angle between the generatrix 
and the axis of the cone 
A = \V/1 + he, coefficient of first quadratic form 
Ni, Ne = elastic normal forces 


= shearing forces 
Z = components of the total load vector 


S,M:, = tangential force and bending moment 
Q: 


The positive directions of forces and displacements are 
shown in Fig. 2. 

Let us make a change of variables which is very convenient 
in the analysis of conical shells (3), using z = ae*', and let us 
rewrite Equations [1] in the following form 


as 
— kN: + + ArX =0 


a 
20 +1rQ: = 0 


+ + ArZ = 0 kM; — ArQ, = 0 


06 


Knowing the internal forces and the moments, we can ob- 
tain the strain components from the following relations 


Translated from Izvestiia Vysshych Uchebnykh Zavedenii 


MVO, Aviationnaya Tekhnika (Bull. Higher Institute of Learning, 4 
Aviation Engng.), no. 2, 1959, pp. 72-82. = Ehe, Ns = Ehe S =0.5Ehw M, = — — {3] 


1 Numbers in parentheses indicate References at end of paper. _ ie 7o 3 Go 
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where 


E = modulus of elasticity of the material 
h = thickness of the shell along the generatrix ros 
h = mh, reduced thickness of the shell along the parall 


The strain components are expressed in terms of the com- 


ponents of the total-displacement vector as follows 


Ar r 00 Adt\r 
1 ov 
1 ow 


‘If we eliminate N2, Q, and Q, from Equations [2] we obtain 
the following system of equations 


a k? Ak 28 as 
(5 +s) m- 
Ar(kZ — X) 
A? 102 a 
aZ 
Ar 42 +7) 
After eliminating S we get ve 
2 
2? 


We shall solve the problem in terms of displacements. 
Using the supplementary condition ¢, = 0, considered in 
Vlasov’s monograph (2), expressions [3 and 4] are rewritten 


‘i Eh du 1 Ou +2 (2) 
Ar ot r 00 
Eh® k? Ou 10 
262" A +5 5(4 


Inserting Equations [7] into Equations [5], we separate 
the variables, eliminate from the system the components of 
the displacement vector u, simplify the coefficients of the 
resultant equation, and get as the final result 


el (1) 


This equation characterizes the elastically stressed and de- 
formed state of a conical shell with a deformable contour. 


If we put in this equation k = 0, and consequently kae** = R, 
we obtain the equation 
proposed by Vlasov (2) for an orthotropic cylindrical shell of 
small length. 
In the foregoing separation of variables, we had in mind the 
Fourier expansions 


u = cos nb X = cos n0 
v = >, sin nO Y = sin nO 
= cos w= cosn6 Z= > cos 


Following Vlasov (1), let us also consider the supplementary 
conditions €, = 0and w = 0. Now the expressions [8 and 4], 
with [9] taken into account, can be rewritten as 

vow 


Eh (= 
Atm \dt 
If Equations [10] are inserted into Equations [6] and the 
— + v=— 


coefficients simplified, we obtain 
d 
at) xX - 


d 


N = cos 
= sin nO 
M: 


Nin 


E-h 


where 


6 12 keg 
ety, = — x > Equation [11] characterizes the elastically stressed and 
deformed state of a conical shell with deformable contour 
An ad the absence of shear deformation. If we put in Equation 
& = 0, and consequently kae** = R, we get the equation 
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© 
4 
= 4 > 
<0 
Pim 
t 
A 4 
| | 
k 
—{2—+k—+k?— Am}Z, | [8] 
2 


proposed by Vlasov (1) for an orthotropic cylindrical shell of Brac | 
medium length 
In the derivation of Equations [8 and 11], starting with the We i : 
condition k < 1, we have simplified the coefficients, which, as : ~f | | 
shown by analysis, is equivalent to rewriting the first equa- Tit] 
tion of [5] in the form s eli 


In the following we shall solve only Equation [11]; the 
homogeneous equation corresponding to it 


— kYy” + = 0 


is solved by successive approximations and written in terms of shat 
Bessel functions of complex argument (Thomson functions) 


= + + CPi + Chu + v,* 


where 


where Dy = - - 3.688p%e—** x 
J 
Dy = kiPy — — — Dy — 3.922pte-** X 
= — 0.751 — 5.075 —— Im Vi 0.898 ¢) Re Io(+/i 0.918 ¢) 


= pe *ReVi ¢) 
= + — analogy with the foregoing, we can, by using the method 
Ba a > _ of successive approximation, obtain an expression for the par- 


= — (Cen + + + + 
where 
= —Rel( Vig) 


= + By + 0.751 + 4.562 Rel.(V/i 0.898 ¢) 
= p? Im 
= + — 4.673p* Im 0.918 


We can next write, on the basis of Equations [3 and 4] 


Ehe-** 
Ni» = k Un 


dt 
Ehke-* 
where 


= —pe-t Re Vi ¢) 
Dx. = + 2Dy + By + 4.097pe- K 
Re V/s 0.898 ¢) 
= Im Vi (Vig) 
=> kt®;; + + 4.284 xX 
Im Vi 11(/i 0.918 ¢) 


Finally, from the homogeneous equation corresponding to 
Equation [12], we get 


+ + + Chu + 8.*) 


“Anta 


a = —2(1 — 


ticular solution of Equation [11], no matter what the right- 
hand side is. 
By way of example we shall consider the case when the 


component of the external-load vector varies as 
Z= Zo- cos nO 


The law of variation of the external load with the angular 
coordinate is shown in Fig. 3. It can be represented ana- 
-lytically by the expression 


where 


= 2(1 — 
= (1+ p)/2 


‘ 


Thus, the right half of Equation [11] will have the form 


at X = Y = Oof : 
Atnk*a? Z[(B + 3)k2 — 


= Eh 


fy In those cases when 6 assumes nonintegral values, the gen- 


eral expression for the particular solution can be written 
where 


1 
= fl? 
(B — — 2m)? — 


We can now supplement the system of Equations [13-16] 
with particular solutions and use this system to determine 
arbitrarily the integration constants for the specified boundary 
conditions. 

As can be seen from the expression for the coefficient ao in 
Equation [11], the foregoing equations are correct only when 
n > 2, but it is easy to check that for the particular values 
n = Oand n = 1 it is possible to obtain from the initial Equa- 
tions [6 and 7] solutions that are in accord with beam theroy. 

As an example, we perform the calculations for the shell 
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re 


shown in Fig. 4. It has the following geometrical character- 
istics 


We shall determine only the second-order stresses (n = 2), 
since a preliminary analysis shows that there is no point in 
going to higher approximations. 

From Equations [11 and 13] we get p* = 
quently 


k = 0.316 
n = 1.0 


To = 7.9 cm 
h = 0.25 cm 
= 160 cm 


1.61, and conse- 


= 2.26¢-*!s 


Let us assume that the shell is cantilever-supported on the 
right end, and has a diaphragm on the left end. The bound- 
ary conditions for this case can be written 


kt = 0 Nix(0) = 0 forkt =2 up(2) = 0 
S.(0) = 0 = 0 
In order to determine these conditions, it is necessary to 
calculate first the functions ®; in terms of certain expressions 
for Un, Un, Nin and S,. However, it must be borne in mind 
that the Thomson functions used here are taken in the form 
tabulated in the Jahnke and Emde tables (5), and differ from 
the functions given by Watson (Bessel Functions) and many 
other authors. 


2415 
‘ | 
i 
- 1000 


\ 
\ 
\ 


o| 034.07 xt 
25 w 100 
| 


mbt the system of Equations [13-17], which we present in 


C2 C; C, Right side 
—3.329 0.109 7.180 —v,* 
—0.991 4.301 —0.218 0.533 —Un* 
—0.790 2.022 1.396 —8.029 —Nin* 
—1.541 6.755 0.966 —7.966 —S,* 
It follows that 
= 5 + Matin®* — MuNin* + MaS,*] 
= [Min * + M22Nin* M2S,*] 
1 
C; = D [—Miwa* + Mostin® — MasNin* + 
1 
= D [Min — + MuNin* — 
where 


A_| 
+1000 maid 
/ 
+500 
01034 | 07 wt 
Fig. 6 


Let us assume that the load under consideration is char- 
acterized by the parameter 8B = 0.5. Then 


= 6.210 Nin*(0) = —0.485 


Un*(2) = 6.450 7 S,*(0) = —1.000 7, 


The arbitrary constants of integration will be determined 


matrix form 


D_ = the principal value of the determinant 
M;; = the minor of the determinant 


In this example we obtain 


q 
= 16.989 C; = —4. 470 7. 


Cz, = 2.424 


q 


We next obtain from formulas [13, 14, and 17] the second- 
— 


ary normal and tangential stresses nd. 


Nu 
7 The final oe results are shown in Figs. 5 and 6. 
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For comparison, the figures also show lines corresponding to 
other values of the parameter 8. 

Let us perform the calculations for the same shell, but 
for different boundary conditions. Let the shell be canti- 
lever mounted on the right end, and let there be on the left 
end either a rigid frame that is not subject to flexure, or else 
a massive ballistic end piece, welded to the shell. The latter 
circumstance is quite important, for although the left side 
may move as a part of cantilever beam, the shell cross section 
cannot distort. In this case the boundary conditions have 
the form 


forkt =0 »v,(0) =0 

u(0) = 0 u(2) = 0 
The functions f, of Equations [13-17] will be calculated (5), 
in analogy with the foregoing procedure, for the argument 


gy = 2.26e-*/2, To determine the arbitrary constant of inte- 
gration, we set up the following matrix 


The right side of this matrix, for a load characterized by a 
_— 8 = 0.5, has the values 


v2*(2) = 5.210 


= 


Npw 


_ values for the arbitrary constants of integration 


aed 41.781 C, = -2.612 


for =2 


C: C; C, Right side 

—4,912 0.790 —2.074 | 

—3.329 0.109 7.180  —v,*(2) 
4.992 —1.541 6.915 —ua*(0) 
4.301 -—0.218 


1000 
— vt 
(AW. 
Fig. 8 


8.425 LG = -0.515 

We then use Equations [13-17] to find the steeees of in- 
terest to us. The results of the calculations are presented 
in Figs. 7 and 8. 

In conclusion let us indicate that the law of distribution of 
norma! stresses along the generatrix (Fig. 5) agrees qualita- 
tively with the results obtained by Kut’inov (4) for a caisson 
with small taper. 

The proposed equations make it possible to determine the 
stresses produced by internal pressurization of the conical 
shell and due to the stiffness of the end frames or supports for 
nonaxially symmetrical loading of the structures. These 
stresses may be quite considerable. For example, the second- 
ary normal stresses for the structure shown in Fig. 5 are 
approximately 30 per cent of the stresses calculated by the 
law of plane sections at » = 0.5 and 8 = 0.5, and 50 per cent 
of the latter at u = Oand 6 = 0.5. 

At larger values of the parameter f, these figures may in- 


- erease, since in this case the load assumes the character of a 
A local load near the point of support. 
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T° FIND a nnaitttalibin relation between the temperature 
T, measured by optical means and the time averaged 
thermodynamic temperature T, it is necessary to know the 
law of variation of the instantaneous value of 7 with time, 
T = f(t). Such a law was derived in a recent experimental in- 
vestigation by N. V. Kokushkin (2). The only temperatures 
observed were in the turbulent combustion zone, temperatures 
T, of the fresh mixture and temperatures 7’, of the combustion 
products. The fraction of intermediate temperatures 7’; is 
negligibly small (7; > T; > 7;). The time averaged tem- 
perature is then 


T = Tip + 


= probability of we a point over 
given time interval 
Pp: = probability of temperatures 72 


The average optical T) is found the 


tions 


T(t) = true temperature alow 
Ce = the constant in Planck’s radiation formula =~ 
radiation wave length 
t = time rr 


where 


Translated from Jzvestiia Akademii Nauk SSSR, Otdelenie 
Tekhnicheskikh Nauk, eg t Avtomatika (Bull. USSR 
Acad. Sci., Div. Tech. ’Sci., Power and Automation), no. 3, 1959, 
pp. 100-103. 

1 Numbers in parentheses indicate References at end of paper. 
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Connection Between the Flame 
Temperature Measured by Optical 


- Means and by Other Methods 


It was found in many investigations that optical methods yield higher values of the temperature 
in the combustion zone of a turbulent flame than other methods. This phenomenon was explained 
by Kadyshevich et al. (1),! who have shown that if the temperature varies with the time at a given 
point, optical methods give values that exceed the time averaged temperature. 

We consider in this paper the use of optical methods (line reversal and infrared pyrometry) to 
measure the temperatures of turbulent flames. A connection is found between the mean thermo- 
dynamic temperature and the temperature obtained by optical data. 

From the values of the combustion temperature as measured by reversal of spectral lines it 
is also possible to calculate the mean square pulsation temperature, caused by insufficient micro- 
mixing of the fuel and air, as frequently occurs in ordinary technical apparatus. 


= 
If relation [1] is satisfied, we obtain from Equation [2] 


1 P2 
exp (c2/AT>) — 1 exp (e2/AT:) —1 exp (e:/AT) 1 


[3] 


Discarding the second term (at 7), the radiation contribu- 
tion is small) we get 
T-T; ) 
T, [4] 


If the value of the average thermodynamic temperature is 
known, the value 7, obtained by a given optical method 
can be calculated (from the value of \). The results of such 
a calculation are shown in Figs. 1 and 2. The calculations 
show that all optical methods give temperature values higher 
than the average thermodynamic temperature, and the excess 
increases with diminishing radiation wave length \ and as ps 
decreases toward the edges of the flame. 

It is more interesting to solve the inverse problem, that 
of determining the average thermodynamic temperature from 
a measured value of 7>. This can be readily done by trans- 


forming Equation [3] 
exp (c:/AT») — 1 


Since 7, is known for a specified value of a of a homogeneous 
mixture, and since 7» is found experimentally, pz can be cal- 
culated. (In all cases, complete chemical combustion, g = 1, 
is assumed.) 

The results of a similar processing of experimental data of 
Kadyshevich et al. (1) are shown in Fig. 3. 

It can be seen that the temperature profiles, calculated from 
values of 7, obtained with an infrared pyrometer, are in good 
agreement with the results obtained by the chemical analysis 
method. Matters are worse if the data on the reversal of the 
spectral lines (Na, Li) are used, owing to the effect of the 
relative error of the measurements. As is known, the total 
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Fig. 1 Comparison of profiles of the mean thermodynamic tem- 
perature with the profiles of mean optical temperatures for vari- 
ous wave lengths }. 1—Profile of thermodynamic tempera- 
ture; 2—Profile of mean optical temperature, obtained from a 
specified value of the thermodynamic temperature, \ = 4.3. 
(measured by infrared radiation); 3—Profile of mean optical 
temperature for \ = 2.74 (measured by infrared radiation); 4— 
Profile of mean optical temperature, \ = 0.6u (reversal method) 


OW 


Fig. 2 Comparison of experimental and theoretical tempera- 
ture profiles in transverse cross section of the turbulent cone of 
a flame; distance y in millimeters from the axis of the cone (a = 
1.2). 1,1‘/—Calculated temperature profile, with the mean ther- 
modynamic temperature taken from the experimental chemical 
analysis data; 2,2’—Calculated temperature profile, with thermo- 
couple readings used for the average thermodynamic temper- 
ature; open triangle—Experimental temperature profile, ob- 
tained by reversal of the spectral lines; open circle—Experi- 
mental data obtained by the method of infrared radiation; 3— 
Experimental data obtained by chemical analysis; 4—Experi- 


error in the determination of pe. 


error in determining the true csineaian by optical methods 
amounts to AJ’ = 50-80 deg. Let us estimate the relative 
We use the transformed 


°s 15 25 A Ap: _ 14,320 50 to 80° 


IN 2 


Equation [5] and the simpler Wien equation 


— 


The relative error is 


(Tr max = 2400 deg abs. for a = 1.) . : 
In the case of the live reversal method 


2400? 0.21 to 0.33, 


06 

for the infrared pyrometer 
= 0.070 to oud 
It is seen thus that a small deviation in the experimenta! 
value of 7», as determined by the reversal of the spectral 


lines, results in a large error in the determination of 
Therefore, even if corrections are included in the calculations 
; ; this method cannot be recommended to determine the profiles 


of the average thermodynamic temperature in the transverse 
_ cross section of the turbulent cone section of the flame. 
Measurements of combustion temperatures by reversal of 


; _ spectral lines in the range 0.8 < a < 1.1 give temperatures 


20 40 60 


Fig. 3. Comparison of experimental and theoretical temperature 
profile in transverse cross section of cone of flame (a = 1.0), 1— 
Profile of mean calculated thermodynamic temperature, in which 
the average optical temperature necessary to determine the 
thermodynamic temperature is obtained from data on the re- 
versal of the spectral lines; 2—Calculated profile of mean ther- 
modynamic temperature, the optical temperature obtained with 
infrared pyrometers; 3—Experimental profile of mean thermo- 
dynamic temperature based on chemical analysis data; 4—Pro- 
file of average optical temperature, experimental, based on data 
obtained by the reversal method; 5—Profile of average optical 
based on optical parameter data 
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Table 1 

Or To AT: ATr To AT, ATr Tor 
1A 2290 2280 2350 0.10 0.09 70 60 2295 0 5 2280 
1.2 2190 2160 2320 0.22 0.13 160 130 2200 35 10 2195 
1.3 2100 2010 2260 0.22 0.15 250 160 2115 35 15 2045 
1.4 2030 1860 2180 0.29 0.16 320 150 2040 55 10 1915 
1.5 1950 1680 2060 0.39 0.12 380 110 1965 140 15 1820 
Notes: spectral-line reversal Ty—by infrared pyrometer 
AT, = Ty — 62, ATr = Ty — 6r—difference in readings between the optical and true temperatures. Oe he! 


that are considerably higher than the possible theoretical 

temperatures at a given average value of a. This result was 

attributed by Kadyshevich et al. (1) to the presence of micro- 

irregularities in the composition of the mixture, i.e., to a 
deviation of the instantaneous values of a from the time 
averaged ones a*. (The fuel benzine was injected through 
a fuel collector made up of 12 nozzles into air preheated to 

200 C at a distance of 4000-6000 mm from the point of 

measurement.) In spite of the fact that, in this case, the 

mixing should have been complete enough for the parameters 

to assume their average values, the actual mixing could not 

be expected to be microscopic. Let us attempt to calculate 

the possible fluctuation in a, or, more correctly, the cor- 

responding temperature fluctuation in these experiments. As 

is obvious from the foregoing, the optical methods measure 

the mean energy of radiation J*, which corresponds to some 

average temperature 6. 

The instantaneous fluctuation of a about the average value 
a* should lead to corresponding fluctuations in the combus- 
tion temperature J 


[7] 
‘as 


@ = average temperature corresponding to a* a 

2 = instantaneous fluctuation of the temperature due to the 
instantaneous change in a 


where 


We assume that the temperature fluctuations are not too 


large, V82/0 < 0.5, and that they are random. Let us find 
the running average value of the radiation energy as a function 
of T. Taking Equation [7] into account, we have 


ov 1 (ov 
IMT) = J*(8) + (=)... 8] 


For J ot we use the Wien formula instead of Equation [2] 


Confining ourselves to three terms in expansion [8], we obtain 


C2 Ce Ce 1 
exp (- x 1.) = exp (- x + ] 
[9] 


We now obtain from [8] the mean square pulsation tem- 
perature V 82/6 = 3*. 

Knowing the value of 7) (measured by reversal of the spec- 
tral lines or by infrared radiation) and the value of 6 we can 
find #*. Two values were computed: 

1 The value of 8, was found from chemical analysis data 
(corresponding to #,*). 

2 @7 and #,* were obtained by thermodynamic calcula- 
tions, assuming complete combustion (g = 1.0). The results 
are listed in Table 1. 

In the right part of the table we give the calculated values 
of T> obtained by infrared methods, assuming that the magni- 
tude of the pulsations of the temperature #* in the stream is 
specified and obtained from the left part of the table. 

The results obtained show that the temperature fluctuations 
due to poor mixing are quite considerable, 15 to 30 per cent, 
and increase with increasing a. 

A measurement of 7, in the combustion product region by 
reversal of the spectral lines and a comparison of this tem- 
perature with the theoretical value (or that obtained by chemi- 
cal analysis) makes it possible to calculate the possible in- 
homogeneities in a and the fluctuations of temperature with 
time. As a result it should also be noted that the use of in- 
frared radiation for the determination of temperature fluctua- 
tions due to inhomogeneities of a with time is not convenient 
in the range of investigated values of a (a = 1.1 to 1.5). 
The difference between JT, — @ is within the limits of the in- 
strumental error of the infrared pyrometer. 
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PTIMAL control systems have recently attracted the 
attention of many investigators (1 to 10),! and principal 
attention has been paid to the analysis of steady-state oscilla- 
tions in the system, since the normal operating condition of 
most optimal contro] systems is the mode of oscillation about 
the position of the extremum. To analyze such a mode of 
steady oscillation in optimal relay systems, use is made of 
expansion of the periodic solution in a Fourier series, in which 
a count is taken either of all terms in the series (4,5) or only of 
one term (8,9). Use can also be made of other methods of 
finding a periodic solution, similar to the investigation of self- 
oscillations in simple relay control systems (11). 

However, as is known from the general premises of control 
theory, in order to obtain a complete and all-embracing esti- 
mate of a control system, it is necessary to consider also their 
transients, and this applies in particular to optimal control 
systems. As will be shown in the following, the investiga- 
tion of transients in such systems can be carried out in many 
cases by the well-known methods of nonlinear control system 
transients as employed in control theory. 

The analysis which follows is limited only to control 
systems in which the extremal characteristic has only one 
maximum. 

1 Let, under certain physical conditions, the extremal 
characteristic shown in Fig. la have the form 1, and let the 
coordinates of the control system be characterized in magni- 
tude by the position of the point of extremum A. Let us 
assume that suddenly the physical conditions of the controlled 
object have changed in such a way, that the extremal char- 
acteristic has shifted and taken position 2 in Fig. la, with a 
maximum point O. Transients will then occur in the sys- 
tem, and will continue until the regulated object arrives at 
the point of new maximum 0, or until steady-state oscillations 
are established about this point. Let us choose the origin of 
the coordinates in the new position of the maximum and 
let us write down a system of differential equations for the 
transients in the optimal system. The initial conditions are 
determined, as usual, from the static equations, assuming 
that the object has occupied the old extremum position A or 
oscillated about this point. We next assume that the change 
in the physical condition was insignificant and that the con- 
trol system has been subjected to ‘‘a sufficiently small dis- 
turbance.” Then the equation of the curve with the maximum 
can be expressed, near the point of maximum, by means of a 
quadratic function. Here as in ordinary control systems, the 
concept “sufficiently small disturbance” remains sufficiently 
broad when investigating the response quality of the tran- 
sients. 

We assume furthermore that the controlled object, which has 
an extremal characteristic, can be represented by means of an 
equivalent circuit (Fig. 1b) in which an inertialess link is con- 
nected in series with a linear link that takes into account the 
inertia of the object. 

2 Let us consider linear systems of extremal regulation 
with a proportional element. The equations of the links of 
the system, as shown in Fig. 1b, have, under the assumptions 
made, the following form. 


Translated from Jzvestita Akademii Nauk SSSR, Otdelenie 
Tekhnicheskikh Nauk, Energetika i Avtomatika (Bull. USSR 
Acad. Sci., Div. Tech. Sci., Power and Automation), 1959, no. 3, 
pp. 32-42. 

1 Numbers in parentheses indicate References at end of paper. 
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The equation of the servo motor is eee aa 
= au (2.1 ] 


The equation of the linear portion of the controlled object 
is 
TE + = Op [2.2] 


The equation of the nonlinear extremal characteristic of the 
controlled object (after the disturbance) is aro 


The equation of the measuring element is iia. 
dp 

where Pee 


u, M, ~,2 = functions of time and represent the coordinates 
that characterize the state of the system 
Q1, = constant coefficients 


The initial conditions are determined from the state of the 
system preceding the disturbance. 

Eliminating the variables u, » and g from Equations [2.1 
to 2.4], we obtain the linear differential equation 


Ts + 20 = 0 [2.5] 


which determines the transients in the system. It is possi- 
ble to continue the analysis of the transients in the extremal 
control system with proportional measuring elements, follow- 
ing the known rules of transient analysis (stability, response 
quality) in ordinary linear control systems. Since the equa- 
tions of the proportional measuring element of type [2,4] 
actually reduce to a linear relation between the coordinates 
z and u, then, if the equations of the servo motor and the ob- 
ject in the control system are linear equations of higher order, 
the optimal control system will nevertheless be reduced to a 
linear differential equation of respectively higher order, and 
can be analyzed in accordance with the general rules. 

It is seen from Equation [2.5] that at different values of the 
amplification coefficients a1, a, a3, the transients have a 
damped character, i.e., the system arrives at the maximum 
point O ast and 2-0. The degree of damping is deter- 
mined here only by the time constant of the object, and the 
frequency also varies as the amplification coefficient is 
changed. 

3 Let us consider next an optimal control system, which 
differs from the foregoing one only in that the proportional 
measuring element is replaced by a relay-type measuring ele- 
ment, the equation of which is of the form . 

Bal] 


The transients in this control system are then described by a 
system of differential equations [2.1 to 2.3] and [3.1]. Elim- 
inating the intermediate variables from this system of equa- 
tions, we obtain a nonlinear differential equation in x 


TE + signz = 0 [3.2] 
The transients which are described by Equation [3.2] are 
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Fig. 1 a—Characteristics of the regulated object; b—structural diagram of optimal control system with direct feedback proportional 
to x; c—structural diagram of optimal control systems with a model of the linear portion of the regulated object: 1—extremal regu- 
lator, 2—actuating motor, 3—linear portion of the regulated object, 4—nonlinear element, 5—linear element, 6—tachometer generator 


best analyzed in phase plane z,v, where v=2. Let us next analyze an optimal control system (Fig. 3a) with 
After integrating Equation [3.2] we obtain the equations of an element of logical action (1,2). 
the phase trajectories (Fig. 2a). Let us consider the case when it is possible to establish 


measurement with respect to z. For simplicity we shall as- 
Jz lees | [3.3] sume that relay 2PR is ideal but relay 1PR is not. Then the 


= + In 
Z = — T(v — + 


The plus sign holds here when z>0, and the minus sign aa Shad ow jo 
when z<0. The symbols and denote the corresponding is 
initial conditions. 

Let us examine an optimal control system which may be i 
described by Equations [2.1 to 2.3] and by the following — 
measuring element equations 


d 


where 2 A is the width of the hysteresis loop of the two-position _ 
polarized relay. Fig. 2b gives the phase trajectories of this _ 
system. 

For the purpose of plotting we used formulas [3.3] and, | 
cording to Equations [3.4], Equation [3.3] holds in the region | 
to the right of the line EFGH when z>0, and to the left when : ‘ 
<0. As can be seen from Fig. 2b, a stable limit cycle is a 
established in the phase plane, and the system returns to it 
ars of the initial conditions. 


= 


‘ake 
= AX 


2 
| Z 


La 


4 4 Fig. 3 a—Structural diagram of optimal control system with a 

“ La logical-action element: 1—logical-action element, 2—actuating 

0 “ motor, 3—linear portion of the regulated object, 4—nonlinear 

element, 5—differentiator, 6—tachometer generator, 1PR, 
Fig. 2 a—Phase trajectories of optimal system with measuring 2PR—two-position polarized relays; b—phase trajectories of 
element described by [3.1], a:a2 = 10, a3 = 0.5, 7 = 1 sec; optimal control ro age with logical action element, aja: = 10, 
b—phase trajectories of optimal system with measuring element a, = 0.5, 7 = 1 sec, A = 0.5 (the dotted lines limit the switching 


[3.4], a2 = 10, ag = 0.5, 7 = 1 sec, A = 0.5 boundary) 


January 1960 


iP 
+ 4 
st ) 
4 
eal, 
2] 
> 
=, 
= | 
e 
/ if 
Zz 
i | 
7 
| 


rf PR—two-position polarized relay 


Fig. 4 a—Principal diagram of regulator with measuring 
element [4.3]; b—principal diagram of regulator with measuring 
element [4.10] ; c—principal diagram of regulator with measuring 
element [4.11]; d—principal diagram of regulator with 
measuring element [4.11], A—amplifier, A M—actuating motor, 


equation for the law of regulation can be written in the an 
u = sign (g + A) sign z [3.5] 

Equation [3.5] can be rewritten, by using [2.3] as 
u = —sign (2a;zv — A) sign v [3.6] 


Since the control coordinate u can assume values of only 
u= +1, then after reduction, the dynamic equation of the sys- 
tem becomes 


TE + Z + Qide signu = 0 [3.7] 


The phase trajectories of this system are given by Equation 
[3.3] and the boundary on the phase plane, corresponding 
to the transition from z>0 to z<0 and vice versa, is deter- 
mined from Equation [3.6]. This boundary comprises the 
two branches of the hyperbola 


v = (A/2asz) [3.8] 


Since the coordinate z is usually not accessible for direct 
measurement, then to obtain the value of z it is necessary to 
use an equivalent circuit for the linear portion of the regulated 
object. The greater the hysteresis loop of the relay 1PR, the 
less stringent the accuracy requirements for this equivalence. 
In individual cases it is merely necessary to introduce time 
delay in relay 2PR. It is frequently better to use a model of 
the linear portion of the controlled object. Fig. 1¢ shows the 
block diagram of a control system with an equivalent circuit 
for the linear portion. 

The value of z is determined from Equations [2.1 and 2.2] 
in operator form 


Consequently, as the time constant of the mndel Tm of the 
link (Fig. 1c) is so chosen that 7,,=7, then 


oa 


F 


Fig. 5 Phase trajectories of optimal system with measuring 
element [4.1], = 10, ag = 0.5, 7 = 1 sec, A = 0.5 


We note that when using an equivalent circuit for the linear 
portion of the control object, it is possible to dispense with one 
of the differentiating elements. 

Using the same method of phase-plane transient analysis, 
it is possible to investigate an optimal control system with a 
memory for the maximum value of the controlled quantity. 
The circuit with capacitor recharging (3) used in such systems 
is essentially none other than the simplest equivalent circuit 
of the controlled object. 

4 Let us consider the improvement in the transient re- 
sponse in optimal control systems, an improvement based on 
introducing into the control law additional terms containing 
the velocity z in some form or another. The advantage of 
such a method is that for normal operation of the control sys- 
tem it is neeessary in general to have the value of z, and there- 
fore the introduction of additional feedback proportional to z 
does not raise any technical difficulties. 

We illustrate this using as an example an optimal control 
system with proportional measuring element, described by 
Equations [2.1 to 2.3] and by the equation of ‘the eed ‘meas- 


uring element 


Eliminating the intermediate variables from these equa- 
tions, we obtain Tye 


r+ (1+ + =0 


By varying the parameter £ it is possible to vary the degree 
of damping. We recall that this could not be done for Equa- 
tion [2.5] for the transients in the system [2.1 to 2.4], to which 
our system reduces when p=0. 

Let us consider an optimal control system which is described 
by Equations [2.1 to 2.3] and the equations of a new meas- 
uring element 


| 
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D 
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The principal diagram of such a diagram is shown in me : a) aN 


4a. Eliminating the intermediate variables from Equations 
[2.1 to 2.3] and [4.3], we get 


+ (1 + Bonas)é + sign (2c + A)=0 [4.4] 


from which we determine the equations of the phase trajec- 
tories 


1 + av/aide 


Here a=1+fmaz; the plus sign holds for 2>(1/2)A/as, 
and the minus sign for z = — (1/2) A/as. 

Fig. 5 shows the phase trajectories for the motion of the 
system based on these formulas, for feedback coefficient values 
B=0.1 and B=0.4. It is seen from this construction that in- 
creasing the coefficient 8 permits not only an increase in the 


damping of the system, but also reduces the amplitude of ia ; 


steady-state oscillations. 

We shall show further, that if the servo motor is relay 
controlled, introduction of feedback proportional to z also 
improves the quality of the operation of the extremal system. 
Let us consider a system described by Equations [2.1 to 2.3] 
and the following equations for the measuring element 


u=sign(u—A) w>0 
u=sign(u +A) w<0 
Here iy 
[4.7] 


Since u can assume the values u= +1, the phase trajec- 
tories for the motion of the system are determined by Equa- 
tion [3.3] with a plus sign when u=1 and a minus sign when 


ma}, 
The switching limits are ~ ss 
4.8] 


u u+A<0 


or, using Equation [4.7], are represented in the zv plane in the 
form of two parallel lines 


—-Bv-A>O [4.9] 

Fig. 6a shows the phase trajectories for 8=0.2 and 8=0.4. 
It follows from this graph that for a relay controlled servo 
motor, feedback proportional to ¢ improves the stability of the 
system and decreases the amplitude of steady-state oscillations, 
At a definite ratio of control system parameters and the value 
of the coefficient £ the transient may become similar to the so- 
called “gliding mode,” in ordinary control systems (11). 

Let us note here a characteristic feature of optimal control 
relay systems. It is known (11) that the introduction of pro- 
portional feeback in ordinary relay systems increases the de- 
gree of its damping. A similar result is obtained by introduc- 
tion of derivative feedback (proportional to %) in an optimal 
control system. 

Let us consider further still another method of improving 
the operating quality of an optimal system. We considered 
before the system [2.1 to 2.4], in which the measuring instru- 
ment performed the function of division. However, realiza- 
tion of such a measuring element is quite difficult. Much 
simpler, in the structural sense, are regulator systems (Figs. 
4b, 4c, 4d) in which the measuring elements have the equa- 
tions 
[4.10] 


[4.11] 


Let us consider first the optimal system of Fig. 4b. Elim- 
inating the intermediate variables from the system of Equa- 


>) 
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Fig. 6 a—Phase trajectories of optimal system with measuring 
element [4.6] (vo = 7, = 0.8), = 10, as = 0.5, 7 = 1 Bec, 
’ : A = 0.5, at 8 = 0.2; b—the same but with6 = 0.4 


tions [2.1 to 2.3] and [4.10], we get _ 7 


TE + & + 2ajQ2a,r2t = 0 [4.12] 


Hence, going over to variables z, v, we obtain after integra- 
tion the phase trajectories (Fig. 7a) 


. The upper sign in Equations [4.12 and 4.13] holds for »>0, 
and the lower one for v<0. 

Comparing the quality of the extremal systems, Fig. 1b 
and Fig. 4b, we conclude that the systems have nearly 
equal dynamic properties, but the system of Fig. 4b is con- 
siderably simpler. 

To improve the dynamic of the system of Fig. 4b, we intro- 
duce feedback proportional to # (Figs. 4c, 4d). The motion in 
the system is consequently described by Equations [2.1 to 2.3] 
and [4.11]. 

Going over to variables z,v we obtain after integration the 
sought phase trajectories (Fig. 7b) 


1+ 
T 


v= — 29) (as [4.14] 

In Equations [4.14], the upper sign holds for »>0 and the 
lower one for v>0. 

Comparison of Figs. 7a and 7b shows how effective the intro- 
duction of additional feedback is. 

5 Let us consider a case in which the external character- 
istic of this system can be approximated, with accuracy suffi- 
cient for practical purposes, by straight line segments 


= sign [5.1] 


a us analyze the transients in a system described by Equa- 
tions [2.1, 2.2, 2.4, and 5.1]. Since 


(dp/dx) = —ay sign x 
u = sign z 


and after eliminating the intermediate variables u and yu, we 
obtain one differential equation 


TE + + sign z = 0 [5.3] 


which indeed describes the transients in the system. 
Equation [5.3] differs from [3.2] only in the magnitude of 

the constant coefficient, and its investigation in the phase plane 

is analogous. Consequently, the system described by Equa- 

tions [2.1, 2.2, 2.4, and 5.1] is analogous to the relay system 

described ‘by Equations [2.1 to 2.3] and [8.1]. 
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~ Let us consider further the transient processes in the sys- 
tem, describable by Equations [2.1,2.2,5.1 and 3.1]. 
Using [5.2], Equation [3.1] can be written 


After eliminating the intermediate variables from Equations 
(2.1, 2.2 and 5.4], we obtain Equation [3.2], which describes 
the transients in this system. Consequently, the optimal 
control system describable by Equations (2.1, 2.2, 5.1 and 
3.1] is analogous to the system describable by Equations [2.1, 
to 2.3 and 3.1]. The transients in the systems differ only in 
the magnitude of the coordinate ¢. 

The addition of a term proportional to the velocity into 
the law of regulation of a system with an object character- 
istic [5.1] also affects favorably the quality of the transients. 
To verify this, it is enough to consider a system with a meas- 
uring element describable by the equation 


dg 


u = sign — Bt = — sign z — Be [5.5] 


Let us examine now an optimal system, describable by 
Equations [2.1, 2.2, 5.1 and 4.10]. 

Differentiating Equation [5.1], we have ¢=as¥ sign z, and 
then Equation [4.1] can be written in the form 


[5.6] 


Eliminating the intermediate variables from the system 
[2.1, 2.2 and 5.6], we get 


TE + sign sign = 0 


u = sign sign 


[5.7] 


which describes the transients in the system. 

Integrating Equation [5.7] in each quadrant of the zv plane, 
we obtain the equations of the phase trajectories in the form 
of straight line segments (Fig. 8a) 


[5.8] 


The quality of the transients in the system improves as 
additional feedback proportional to z is introduced. In this 
case the system is described by Equations [2.1, 2.2, 5.1 and 
4.11], whereas the phase trajectories are straight lines (Fig. 
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Fig. 7 a—Phase trajectories of optimal system with measuring 

element [4.10], a:a2 = 10, a: = 0.5, 7 = 1 sec, A = 0.5; b— 

phase trajectories of optimal system with measuring element 
[4.11], = 10, a3 = 0.5, = Lsec,A = 0.5 


r(v — %) + + mala, + [5.9] 


The upper sign in Equations [5.8 and 5.9] corresponds to 
one of the combinations of inequalities of the same sign (v> 
0, z>0) or (v<0, z<0); the lower one corresponds to one 
of the inequality combinations of opposite signs (v>0, z<0 or 
v<0, z>0). 

In practice, the case where the extremal characteristic is not 
symmetrical and can be approximated in the form of straight 


"4 


Fig. 8 a—Phase trajectories of optimal system with measuring 


element [4.10] with the characteristics of the regulated object 

being approximated by Equation [5.1], azz = 10, ag = 0.75, 

t = 1 sec; b—phase trajectories of optimal system with measur- 

ing element [4.1] with the characteristics of the regulated object 

being approximated by Equation [5.1], aia: = 10, a3 = 0.75, 
= 1 sec, 6 = 0.4 
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line segments having different slopes 


sidered, additional feedback proportional to z improves con- 
siderably the dynamic properties of the system. Actually, 


> ae gg = —Auzr [5.10] solving the system of Equations [2.1 to 2.3 and 6.2] together 
ae with the equation of the new measuring element 


when z<0, is much more frequent than an extremal 


istic approximated by a parabola. In this case the investiga- 
tion of the transient in the optimal control systems can also 
be carried out in the phase plane, as was done by Ostrovskii et 
al. (10). 

6 In the papers of Kulebakin (13), and later by others, it 
was shown that it is possible to improve the dynamics of 
ordinary servo systems by intoducing error derivatives. This 
method can be extended also to optimal control systems. 
This results not only in improved system quality, but consider- — 
able simplification of the investigation of the optimal systems. — 


Let us consider systems in which the object of regulation © 


can be represented by an inertial-less nonlinear link in series 
with an inertial linear link. 

Let the optimal control system be described by Equa- 
tion [2.1] for the actuating motor, by Equation [2.3] for the 
object, by 


9 ny + y= asp [6.2] 


and by the followi ing equation for the measuring element 


u= a! (6.3) 
Solving Equations [2.1, 2.3 and 6.1 to 6.3] simultaneously, 
we get 


a) 
Qné+ + [64] 
This equation differs from Equation [2.5] only i in the co- 
efficient preceding the highest order derivative. 
Introducing a second derivative into the measuring element 
of this system, its equation becomes 
% 


After eliminating the intermediate variables from. Equa- 
tions [2.1, 6.1, 2.3, 6.2 and 6.5] we obtain the first-order equa- 
tion 


165) 


Thus, the order of the equation has been reduced, and the 
quality of the system has improved. 

Let us consider now an optimal system in which the object 
of regulation, jointly with the inertial measuring element, can 
be represented by a series-circuit comprising a linear, nonlinear 
and a linear link, and describable by a system of equations 
[2.1 to 2.3, 6.2 and 6.5]. 

Solving these equations simultaneously, we obtain a linear 
differential equation for the dynamics of the system 


Z+ = 0 


TE + = 0 [6.7] 


analogous to Equations [2.5 and 6.4]. The introduction of 
the second derivative into the measuring element [6.5] has 
made it possible to reduce the investigation to that of a 
linear equation. Were we to consider the same system but 
with a measuring element [2.4], the dynamic equation of the 
system would be a third-order nonlinear differential equation, 
which is, in general, difficult to analyze. 

It is easy to show that for the optimal system just con- 
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we obtain a silnonas equation of the system 
TE + (1 + + = 0 [6.9] 


in which it is possible to change the damping of the system 
as desired by varing the coefficient 8. A similar optimal sys- 
tem, but with a measuring element of the relay type, describ- 
able by the equation 


= 
& = a sign (6.10] 

is reduced to a dynamic system equation ; 
Te [6.11] 


analogous to Equation [3.2], the phase-plane solution of which 
was obtained above. It is easy to show also that introduction 
of an additional term — fz into the control Equation [6.10] 
makes it possible, all other conditions remaining equal, to 
change the damping of this system. 

Systems in which a linear portion having a high order is 
connected past the nonlinear element are describable by com- 
plex nonlinear differential equations. It is easy to show that 
by introducing into the measuring element of such a system, 
in a suitable manner, corresponding higher derivatives of the 
output coordinate of the object of regulation, will make 
it possible to reduce respectively the order of the overall 
dynamic equation of the system and thereby simplify the in- 
vestigation of the system and improve the transient response. 

7 Let us consider a system of extremal control with con- 
stant time delay, in which the linear portion of the object is 
described by the equation 


+ x(t) = — Ts) [7.1] 


and the remaining elements of the system are described by 
Equations [2.1, 2.3 and 2.4]. Eliminating the intermediate 
variables, the equation of the system becomes 


T(#)(t) + 2#(t) + — = 0 [7.2] 


Thus, an investigation of the dynamics of such a control sys- 
tem can also be performed by the known methods of general 
control theory. The improvement in the transient response 
in this system is attained by introducing additional velocity 
feedback proportional to %. In this case the transients are 
determined by the dynamic equation 


TH(t) + Z(t) + Bayart(t — 73) + — = 0 [7.3] 


and a comparison of this equation with Equation [7.2] can 
demonstrate the suggested improvement in system response. 
An investigation of a similar system, but having a relay- 
type measuring element, such as describable by Equations 
[2.1, 6.1, 2.3, and 3.4] can be carried out in phase plane. The 
system is described by the equations 


+ + aia, = 0 
a(t Ts) >- A 


T#(t) + 2(t) — = 0 
a(t = Ts) <A 


the integration of which yields Equations [3.3], and the switch- 
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Fig. 9 Phase trajectories of optimal system with delay in the 


. object and with measuring element [3.4], a:a2 = 10, 


a; = 0.5, 7 = 1 sec, 73 = 0.25 sec, A = 0.5 


ing limits are of the following form in this case 


AF +2 +1 (exp 2-1). 
T T T 
[7.6] 


The phase trajectories of such a system are shown in Fig. 9. 
The introduction of the derivative z into the measuring ele- 
ments also improve the response of the system. 

If the element with the constant lag is connected in the 
control system past the nonlinear element, something that can 
occur physically if the lag is built into the measuring element, 
then the optimal system is describable, for example, by 
Equations [2.1 to 2.3] together with (t)=¢(t—rs). 

If a measuring element of type [2.4] is used the investiga- 
tion of the system becomes quite (On other 


or 
ing 


¥ 


8 


hand, if we construct a measuring element with an equation 
of the form 


[7.7] 


where z,=2z(t—r3), then after eliminating the intermediate 
variables, the dynamic equation of the system has the form 
[7.2]. To obtain the function x it is enough to introduce into 
the model of the linear portion of the regulated object (element 
6, Fig. 1c) an element having a constant lag equal to 7s. 


u= 
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Schuler’s period (84.4 min). 


not less than two channels. 


I Methods of Damping Inertial Systems 


ET US assume that an object moving at a constant alti- 
tude above the Earth (which is assumed to be spherical) 


carries a gyro platform which is free in azimuth. Let, at the 
initial instant of time, the plane of the gyro platform be 


form two accelerometers with mutually perpendicular axes, 
lying in the plane of the platform (Fig. 1). We orient the z 
and y axes of the coordinate system along the axis of the 
accelerometer, and the z axis upward along the local vertical. — 
We denote by a and # small deflections of the gyro vertical 
from the local vertical, along the z and y axes, respectively, 
and by V, and V, the components of the inertial velocity of 


perpendicular to the local vertical. We place on the gyro plat- a 


the object. If the gyro platform is perpendicular to the local — 


vertical, then the accelerometers mounted on it will measure ~ 
the components of absolute acceleration, double integration of — 
which yields the displacement components. To maintain the 


Translated from Izvestiia Akademii Nauk SSSR, Otdelenie 
Tekhnicheskikh Nauk, Energetika i Avtomatika (Bull. USSR 
—_ = Div. Tech. Sci., Power and Automation), no. 3, 1959, 
pp. 11-18. 
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As is known (1,2)! inertial systems based on the use of the gyroscopic model of the pendulum of 
Schuler have the following principal shortcomings: 1 If the compensation conditions are exactly 
satisfied, the system becomes insensitive to disturbances caused by maneuvering of the object 
carrying the inertial system, but in the case of non-zero initial conditions undamped oscillations 
occur in the system (the system is at the stability limit); 2 in the absence of external correction 
for the drift of the gyroscopes, large errors accumulate in the system with time. 
tions of the vertical of an inertial system could be reduced to zero by introducing damping in the 
system, i.e., by introducing into the equations of motion such terms for which the motion changes 
from periodic into damped. However, the well-known error rate method of introducing damping 
by using internal signals leads to violation of the compensation conditions and does not make it 
possible to obtain rapid damping of the processes. The damping time is comparable here with 


The initial devia- 


There is nothing unique about these features of inertial systems, for they are characteristic of all 
single-loop closed systems. It is known that in single-looped systems it is impossible to satisfy the 
conditions of insensitivity to disturbance (compensation) without violating the stability conditions. 
To satisfy this insensitivity condition at a given point of the system with respect to definite dis- 
turbances it is necessary that the information concerning this disturbance arrive at the point over 
If the signals—the carriers of information concerning the disturb- 
ances—arriving over one channel are in phase opposition with respect to the signals arriving from 
another channel, compensation of signals can occur, and the system becomes insensitive to dis- 
turbances. Consequently acceleration insensitive systems, having acceptable dynamic properties, 
should by their very nature be multilooped systems. The additional loops over which the compen- 
sation signals are transmitted can be either closed or opened. 

This paper contains the theory and principles of construction of damped inertial systems with 
arbitrary period, insensitive to acceleration with respect to the maneuvering of the object. A 
theorem is formulated here, which is a generalization of the corresponding theorem by Schuler. 


gyro vertical in a position perpendicular to the local vertical, 
it is necessary to apply to it correcting signals, proportional 
to the integral of acceleration signals. Thus, the equations 
of motion of the gyro platform will be 7 


R 


where 
ki, ke = gains of the integrating elements (Fig. 2) : 
7 R = Earth’s radius 
g = acceleration due to the central forces acting on the 


gyro platform 


Since Equations [1.1] have an identical structure, we shal 
consider only one of these. 

In the derivation of Equations [1.1] it was assumed that 
the gyro platform rotates about a vertical axis with a velocity 
(wg + i) sin ¢ where wg is the angular speed of the Earth’s 
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‘rotation, the local latitude and the rate of change of 


longitude. 
—s:, If we differentiate the first equation of system [1.1] with 
a ; respect to time and rearrange terms, we obtain 
1 
(p? + ghikeja = (iat Vs [1.2] 
4 —Gyro Platform It is shown in the block diagram of the system (Fig. 2) that 


the accelerometer measures the acceleration Vz — ga, which 
after a single integration yields a signal proportional to the 
29 velocity. This signal is fed, on the one hand, to the correcting 
mechanism of the gyro platform, which has a transfer function 

7 k./p, and on the other hand to an integrator with a transfer 

. uv oft function Rk./p, to obtain the component of the elapsed dis- 


tance. 
If we satisfy in Equation [1.2] the compensation condition 


aout Alen aes If condition [1.3] is satisfied, the inertial system is undis- 
turbed by the acceleration of ‘the object. However, under 
initial conditions (a ~ 0, 0), undamped 
Pik. “ : oscillations arise in the system with a frequency % and with 
an amplitude determined by the initial conditions. If the 
g instrument errors are taken into account, the angle a is 
subject to additional errors. 
Sep tei =f If the condition a = 0 is satisfied, the first integrals of the 
7 Rp? acceleration components represent, with allowance for the 


Earth’s peripheral velocity, the components of the ground 
velocity of the object, whereas the second integrals represent 
the components of the elapsed path. Consequently, it is 
_ possible to determine the velocity vector, the position co- 
ordinates of the objects and other navigational elements. 
To eliminate the undamped oscillation of the inertial sys- 
tem it is necessary to incorporate in it devices that introduce 
damping of the oscillations. Fig. 3 shows a block diagram of 
an inertial system in which the damping is produced by in- 
cluding the first (and perhaps also the second) integrator in a 
positive (dotted) or negative (solid) feedback loop. 
_ The equation of motion of the inertial system with feedback 
damping has the form 


(pt + + Meta = (24 = [1.5] 


If the damping is produced by positive feedback, the equa- 
tion of motion will be 


(p? + 2dQop + = (2d% = [1.6] 


follows from Equations [1.5 and 1.6] that damping by 


ul] Rk, | Se ‘eral inclusion of the first or second integrator in some sort of feed- 
P P h ee back loop destroys the insensitivity of the system [1.3]. It is 
~ 2 oo ee possible to prove the following premise, which is valid in gen- 
ee ® eral: Damping of an inertial system by means of any internal 

; coupling destroys the acceleration insensitivity conditions. 
Sods ae] we It is interesting to note that if the inertial system in which 
g of the damping is obtained by internal feedback is incorporated 
itt mili Aina att in a closed control loop and made to serve as a position-coor- 
dinate transmitter, its damping disappears. In fact, assume 
Rp? 4 that the coordinate n is fed to an autopilot (Fig. 4). The 
2, aL transfer functions of the system in response to a gust of wind, 
iota nas aie = U/,V, and in response to a signal of a change in setting 
94 ARS JournaL SUPPLEMENT 
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of 


It is easy to see that 
A = + — — < O > 


assume the form 
aop(p* + 2dQep + 
PX p* + + M*) + (ap + cep? + M*) 


(cp + cap? + 2dQep + 
No + + + (ap + cp? + %?) 


where the binomial cp + c, characterizes the law of control of 
the autopilot and the dynamics of acceleration in response to 
an error signal applied to the control system. 

Since the coefficients of the polynomial in the denominators 
of expressions [1.7 and 1.8] are positive, the stability of the 
system can be estimated from the sign of the determinant _ 


[1.7] 


|< 


[1.8] 


1 + Ce [1.9] 
0 + 2d% 


[1.10] 


Consequently, if the inertial system is included in a closed 
control loop as a coordinate reference, the damping introduced 
in it by internal feedback vanishes. The closed system thus 
becomes unstable. 

To estimate the character of the divergent transient in the 
closed system, let us consider in greater detail the denomina- 
tor of the transfer functions [1.7 and 1.8], which can be 
approximately factored into 


3 
pt 0.) (p? + + es) [1.11] 


It follows from this that the system is subject to two mo- 
tions, a slow one with the frequency of the inertial system, 
and a fast one with a frequency of the center of mass of the 
aircraft. These frequencies differ by two orders of magnitude. 
Furthermore, the inertial system is divergent, while the de- 
gree of buildup of the oscillations 2d{%?/cz. does not exceed 
several tenths of a per cent during each period (of 84.4 
min). 

Thus, damping introduced into an inertial system by 
means of internal feedback, leads to a violation of the accelera- 
tion insensitivity condition if the inertial system is used as 
an indicator, and disappears entirely if the inertial system is 
used as a control system. It follows therefore that other 
damping methods must be employed. 


2 Use of External Information for Damping and 
Altering the Period 


Let us assume that the inertial system placed on the object 
receives extraneous? discrete (i.e., intermittent in time) in- 
formation (Fig. 5) which characterizes the actual ground 
position of the object, its velocity, its acceleration, etc. Over 
short intervals of time, the extraneous information, which may 
be received from a radio system or from distance-computing 
systems, may have a high degree of reliability (accuracy). 
This information can be used to correct the errors of the 
inertial system, to introduce damping and to change the period 
of the system. Let the information be introduced into the 
system via the input of the first integrator (z,), the input of 
the second integrator (z2) and the input of the gyro platform 
(zs). It may also pass first through special filters. 

The equation of motion of the system (Fig. 5), with ac- 
count of the external information, can be represented in the 


2 Extraneous information is meant to be information on the 
velocity and coordinates of the position of the object, obtained 
by means of any system that is not inertial. 
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following form 


(p? + gkike)a = kykez; + + {2.1] 


It was assumed in the derivation of this equation that the 
compensation condition [1.3] was satisfied. 

Let us choose the externally derived sources of information 
4, 2, and zs such that the vertical’s own motion is damped, 
and the period of the system is arbitrary and differs from the 
Schuler period, while the right half of Equation [2.1] differs 
little from zero. For this purpose we put, for example 


a = —k(a + Ag) — + Ag) a=0 =0 [2.2] 


In expressions [2.2] the angle a and the angle of velocity & 
are determined from 


R 
where ar 
Vz, 8 = velocity and the coordinate of the object measured 
by means of the inertial system 
V0, 80 = velocity and coordinate of the object measured by 


means, say, of radio devices (external informa- 
tion) 

A., Ag = the errors in the determination of the angle and 
angular velocity 


It is assumed that the external information on the coordinate 
and velocity of the object has higher accuracy than the in- 
formation produced by the inertial system. 

If we insert Equations [2.2] into [2.1], we obtain 


(p* + 2dQoxp + Qorx*)a = FM%xX? — F [2.4] 


Depending on the sign of the coefficient k, it is possible 
either to increase (k < 0) or decrease (k > 0) the period of the 
inertial system as compared with the Schuler period, without 
violating the compensation conditions [1.3] (with accuracy 
to within the errors of the system, due to inaccuracy in the 
external information). Since Equation [2.4] contains a damp- 
ing term, the vertical’s own motion will be damped. 

Thus, we can formulate the following theorem: 

Theorem—To construct a damped inertial system with 
arbitrary period, invariant with respect to maneuvering of the 
object, it is necessary to introduce into it extraneous informa- 
tion concerning the velocity and coordinates of the object, 
and to determine the angles and angular velocities of the 
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vertical on the basis of this information and the corresponding 
information produced by the inertial system. 

The external information need not be continuous, but can 
be applied only at brief time intervals, for example during the 
initial display (to eliminate the initial errors), to eliminate the 
accumulated errors after a prolonged operating time, etc. 
The time during which the external information is supplied 
should be as short as possible; during this time the transients 
due to elimination of the errors should essentially be ter- 
minated. 

In order for the transients in the system to be terminated 
within the shortest possible time, the period of the system 
should be considerably less than that of Schuler’s period (x? > 
1). However, a considerable reduction in the period makes 
the inertial system sensitive to noise due to the external in- 
formation. The choice of an optimum coefficient of period 
reduction d is one of the principal problems of the theory of 
inertial systems with arbitrary period. The choice of the 
magnitude of the decrement d is dictated by the requirement 
that the processes must be damped rapidly, thus calling for, 
say, d = 0.7 to 1.0. ‘os ee 


Arbitrary Period 


$3 Analysis of Errors of Inertial System With ip 


The errors A, and A, in the extraneous information, used 
for the correction, give rise to errors in the inertial system. 
It is obvious that if A, = Ag = 0, the intrinsic errors of this 
system will vanish after the transients have died down. 
However, it is impossible to determine by any method what- 
ever the velocity and position of an object without an error, 
and therefore we shall consider the characteristics of inertial 
systems in the presence of the foregoing errors. 

We consider two typical cases: 

1 The errors A, and Ag are independent of the time and 
are constant. 

2 The errors are random functions of time. 

In the case of constant A, and Ag, the steady-state error 
a? of the inertial vertical will be a 
1 2d 

Aa ¥ Ox Aa [3.1] 

It follows from this that at large numerical values of the 
period reduction coefficient x(x >> 1), the major source of 
steady-state error of the system will be the errors in the co- 
ordinate information, whereas at small numerical values of 
x the decisive errors will be those in object velocity data. 

If the inertial system is corrected without changing (reduc- 
ing) the period, the errors in the object velocity information 
may prove to be inadmissibly large. This results from the 
following considerations. Modern navigational systems 
make it possible to determine the object position over ranges 
of several hundreds of kilometers with an accuracy up to 
Asz = 1 km, while the speed can be determined, with the aid 
of a Doppler system, with accuracy to 0.3-0.4 percent. Con- 
sequently 


= 


Aa _ 1 0.6’ Ae 0.004 V0 
R 6 X 10° R 
di: ., taking the object velocity to be Vo = 1 km per sec - 
0.004 x 1 angular min a. 
| ‘ 
Aa = 0:23 X 10 ( ) 


Knowing these errors in the external information, we can 
find the errors in the initial system under various reductions 
in the period. For x = 1 we have 


y % 6 shows plots of the functions Asg. 


- are random functions of time. 


When the period is reduced by x = 10, we find 
= 0.6’ ¥ 0.37’ 


i.e., the error does not exceed 1 angular min. 

Let us estimate the system errors if the errors A, and A; 
After estimating these errors 
in terms of their spectral densities As,(w) and Asg(w), we can 
write for the mean-squared error of the vertical 


= 
8 


The spectral densities of the errors As,(w) and Asg(w) are 
determined essentially by the properties of the systems used 
to measure the quantities s.9 and Vzo, and by the vibration 
of the objects on which these systems are mounted. The in- 
fluence of the object vibrations on the magnitude of the error 
is particularly substantial when the results of the measure- 
ment depend on the position of the object relative to the local 
vertical. In such cases the oscillations, which are random 
functions of time, determine to a considerable extent the 
spectral characteristics of the errors Ad, and Asg. A Doppler 
system for measuring ground speed is indeed such a system 
in which the errors are determined by the vibration of the 
object. When the position coordinates of the object are 
measured with the aid of radio-navigation systems, the effect 
of vibrations is less noticeable. Taking the foregoing into 


- account, we can assume for the spectral densities Asg(w) and 
As,(w) the following expressions (3,4) 


+ w?*) 
As(w) = + w)? — = 6,2 + a*) [3.6] 
where 


= resonant frequency 
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Inserting the expressions [3.3-3.6] into [3.2] we obtain 
after integrating and transforming seal 
where é,¢ and é,, are the mean-squared values of the errors 
of the inertial system due to inexact object coordinate and 
velocity information. With this we have 


x? — 1)X2dMx + fr) 


fea! = + + 
= [3.9] 
where 
A = — + (2h — — 
Ag = — Doty!) + (2h — — 


A, = — Doty!) + (2h — A*) + — 
Ap = 24 — + 2X 2h — 
n = Vogt—1—2dV a? — 1+ 


2h = 2Q%x%2d? — 1) 
= 222 — 4a? 


Let us analyze the errors é, and é,¢. It follows from Equa- 
tion [3.8] that when x = 1, i.e., in the case of a system with 
a Schuler period, the position error é~,, vanishes. As x in- 
creases, the error é,q increases, and as the period of inertial 
system is further reduced, it tends to a maximum é,¢ max? = 

The velocity error also depends on x. If we assume that 
d = 1 and 8; < Q, then as x ranges from 1 to infinity the 


error 
8 


= 8 


Q 


starts increasing, and reaches a maximum 


.2 
= on for x = 1 
o 
2 
Ca 
= 12 for x = 
Upon further increase in the coefficient x the error begins 
to decrease, reaching zero at x > ©. Such a character of the 
dependence of the error é,4 0n x is explained by the fact that 
as the frequency of the inertial system increases, the ampli- 
tudes of the frequency characteristic decrease in inverse pro- 
portion to x?, with a total decrease in the spectral density at 
frequencies above resonance. One must not conclude from 
this that it is necessary to increase the frequency of the inertial 
system to values exceeding the resonant frequency w. When 
the frequency of the inertial system is so increased, additional 
errors due to angular vibrations of the airplane arise. 
The choice of the coefficient x can be based on the require- 
ment x < wo. If wo is equal to 0.02 sec, then x is less 


than 16.3. Apparently, a reduction by a factor of 10 is 
acceptable. 

If the damping parameter £; is greater than Q, the function 
2,4 does not have a maximum with respect to x and decreases 
monotonically with increasing x. 


4 Influence of the Object Speed on Operation of 
the Inertial System 

It has been assumed up to now that the acceleration g due 
to the central forces is constant and independent of the 
velocity V of the object, something that is true at low veloci- 
ties (V < 3c, where c is the velocity of sound). For large 
velocities the acceleration g, to which the elements in the 
object are subjected, will be 


R 
where 


go = acceleration due to gravity 
V absolute velocity of the object 


It follows therefore that for a velocity 
V 


the period of the inertial system tends to infinity. Conse- 
quently, the inertial system cannot function on an object 
that moves with the satelloidal velocity. However, if the 
initial errors are eliminated by using external reference infor- 
mation to reduce the period and the damping, then a large 
period (infinite in the limit) after cessation of the correction 
_ may prove to be a favorable property of the inertial system. 
Thus, to construct damped inertial systems with arbitrary 
period it is necessary to have external information, which 
ean be received by the system either discretely or con- 
tinuously. In the case of discrete information, the pauses 
_ may be arbitrary. One can also conceive of inertial systems 
in which the extraneous information is used only during the 
initial instant of operation for initial adjustment (to eliminate 
initial errors). Such systems then continue to function as 


Schuler systems. 
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: Properties of Air in in 
“the Temperature Interval From 1000 


12,000 K and the Pressure 
Intervals From 0.001 to 1000 atm 


ODERN theory of jet propulsion deals with high tine 
peratures arising when bodies move through the at- 
mosphere at supersonic speeds. In this connection it be- 
comes necessary to study the properties of air under such 
conditions. There are at present no fully developed experi- 
mental methods for direct investigations of the properties of 
gases at high temperatures. Theoretical calculations there- 
fore become particularly important. Methods of statistical 
physics, in its quantum-mechanical formulation, along with 
spectroscopic data, serve as a reliable base for a theoretical 
calculation of the thermodynamic properties of gas systems 
and especially of mixtures of gases capable of reacting 
chemically with each other. An example of such a mixture 
is air: In addition to the dissociation of 0. and Nz and the 
formation of NO, ionization of both the molecular and the 
atomic components occur in air at high temperatures. 

The thermodynamic functions of air are determined from 
the thermodynamic parameters of its components and from 
its composition, which varies with the temperature and 
pressure. Therefore the problem breaks up into two stages: 
Calculation of the thermodynamic parameters of the “pure” 
components, and calculation of the composition of the air 
and ‘its thermodynamic functions. 

The calculation of the thermodynamic parameters of 
“pure” components reduces to the calculation of statistical 
sums. For example, the free energy of one mole of a one- 
component ideal gas is expressed in terms of the statistical sum 


Q'of the internal degrees of freedom of the gas in the following 7 


manner 
F = —NkT ln ( ) Q 
v here = 


statistical weights of the levels 


_ N = number of molecules in gm-mole 
k = Boltzmann’s constant 
‘on = Planck’s constant ay 
T = absolute temperature 
V = volume of gas 


@ = base of natural logarithms 


To determine the composition of the air it is necessary to 
solve a system of nonlinear algebraic equations, which include 
the equations of the law of effective masses for each of the 
possible reactions in the air and for each ionization process, 


Translated from ‘Physical Gasdynamics,” USSR Acad. Sci., 
pp. 3-38. 


the equation for Dalton’s law, the equation of material | 
balance, and the law of charge conservation. 

To calculate the thermodynamic functions of a mixture of 
reacting gases it is necessary to solve two systems of linear 
equations, 


1 Methods of Calculating the Statistical Sums 
for Atoms, Molecules and Their Ions 


The statistical sums of atoms are the simplest in form. 
Henceforth we shall understand “statistical sums” to mean the 
statistical sums of the internal degrees of freedom. 

In caiculating the statistical sums, the energy of the un- 
excited (ground) state of the atom was assumed to be zero; 
i.e., we calculated statistical sums of the form a: 

€1 €2 
Q = auc(go + ge AT 4 

For most atomic systems, several of the first excited levels 
are located 20,000 cm above the ground state, and many of 
the excited levels are located above this small group of levels, 
100,000 cm from the ground state. At low temperatures, 
up to 6000-8000 K, it is possible to neglect the influence of the 
large number of levels located 100,000 cm! above the ground 
level without introducing a noticeable error in the value of 
this statistical sum. However, at higher temperatures (from 
10,000 K and up) such neglect introduces a noticeable error in 
the values of the statistical sums and their first and second 
derivatives with respect to the temperature. Thus, for ex- 
ample, for nitrogen at 7 = 12,000 K, neglecting the influence 
of all the levels located above 30,000 cm—', we would introduce 
an error of 0.6 per cent in the value of the statistical sum, an 
error of 1.2 per cent in the value of the first derivative with 
respect to the temperature, and an error of 18 per cent in the 
second derivative with respect to the temperature. 

In calculating the statistical sums of atoms and their ions, 
we took into account all the known levels, but instead of 
using the exact values of e; of the excited levels, we used values 
rounded off to three significant figures. The error introduced 
thereby in Q is estimated in the following manner. We repre- 
sent ¢; in the form 


is 
= €> + Ae 


where 


The statistical sum becomes 
Q= Le 

p 
Expanding exp (—Ae,;®/k7T) and retaining the first two 
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terms of the expansion, we obtain 


To calculate the thermodynamic functions it is convenient 
to use the quantities Q and 
dQ 


The statistical sums for diatomic molecules are of the 
form 


_ v=vmax ep(l,v) J=J max _ er(lv.J) 
Q=gaue ge *T Die (QJ +1)e 
l v=0 J=J min 
[4] 


where ¢,, €, and e; are the rotational, vibrational and elec- 
tronic energies of the molecule. In the simplest case, when 
the resultant electron spin and orbital-momentum projections 
on the axis joining the nuclei vanish, the rotational energy is 
given by 


= +1) — + + 
+1) + HIT + 1+... [5] 


he 


= 2d,T? + 6fv,T? + 12h,T* + 20m,7'§ — 


J is the rotational quantum number, which assumes all per 
missible values up to Jmax(v, 1). Jmsx(v, l) is the maximum 
rotational quantum number at which the molecule is stable. 

B., D., Fe, He, Oe, Be; Ye and 6, are constants determined 
from the analysis of the spectra. The vibrational energy is 
expressed in the form of the following series 


Ex = Wo — + woyov? + Woz! +... [7] 


where is the vibrational quantum number, and woyo 
and wp2 are spectroscopic constants. 
In the calculation of the rotational portion of the statistical 


sum 
J=J max 


(8) 


the summation over the rotational quantum numbers is re- 
placed by integation. As a result of this replacement, Q, is. 
expressed as a series in powers of 7’ with diminishing coeffi- 
cients. Neglect of the sixth and all other terms of this series 
does not introduce errors greater than 0.001 per cent in Q,, 
0.008 per cent in Q,’/7,, and 0.3 per cent Q,”/T. 

Estimates show that at temperatures above 8000 K a notice- 
able contribution is made to the statistical sums of the di- 
atomic molecules O2, Nz, NO, O2+, Nz+ and NO+ by the energy 
levels located above the limit of dissociation of the electron 
ground state of the molecule. In this connection, we took into 
account in Q, the finite number of rotational levels. Integra- 
tion of the rotational portion of the statistical sum was carried 
out in the limits from Jmin to Jmax. Up to temperatures of 
6000-8000 K, we can integrate from Jmin to ©. For molecules” 
made up of identical nuclei Nz, 02+, N2*), we took into 
count the factor one half in front of the rotational component of | 
the statistical sum. The lower electronic states of the mole- 
cules O2, NO, N2+ and O2+ have electronic terms that are differ- 


ent from ‘>. Therefore, because of the interaction of the 
spin with K or with the axis, splitting of the rotational levels 
takes place. Estimates have shown that this splitting, in the | 
temperature interval from 1000 to 12,000 K, can be ac- 
counted for by means of a statistical weight for Q,, calculated 
under the assumption that there is no splitting, as was indeed 
done in the calculation of Q, for all the electronic states. 
The final formulas for Q,, Q,’ and Q,” are 


beT + 2dyT? + 3f,T? + 4hyT* + — + 1)leT’ + (pe + 2)peT? + 


(pe + + + 4)teT* + (pe + 


he 
+ (ps + 2p, + 2)peT? + 


(Pe? + + + + + 12)teT* + (oe? + + 


where 
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h, = + 12054 + 1050 
D, D, k D4 
= {120 + 1680 ™, = 1680 
he 
= +1) m= 
The vibrational component of the statistical sum 


together with its derivatives, was calculated by direct sum- 
mation over the entire interval of temperatures from 1000 to 
12,000 K over all possible v from 0 to vmx. For the vibra- 
tional component of the statistical sum and its derivative we 
wed the following formulas 


i % In the calculation of the electron portion of the statistical 
sum 


1 
we took into account all the known levels located 100,000 
cm~' above the ground state of the molecule. The equations 
for the electronic component of the statistical sum are 


2 Spectroscopic Constants Used in Calculation 
of Statistical Sums 


? The spectroscopic constants used for the calculation of the 
atomic and molecular statistical sums are listed in Tables 
l1to7. The sources used in the compilation of this table were 
essentially in references (1 to 3).!_ We also used the results of 
later works, which were not included in these collections. 
The source reference from which the corresponding constants 
are taken is indicated in each case. The quantities in par- 
entheses were calculated either by solving the wave equation 
for the given case, or from empirical or semi-empirical form- 
ulas. The tables list all the atomic and molecular electronic 
levels taking into account the calculation of the corresponding 
statistical sums. 

Four, three or two terms of the series of Equation [7] are 
usually given for spectroscopically investigated electronic 
states of molecules. The possibility of determining the 


1 Numbers in parentheses indicate References at end of paper. 


terms of this series depends on the number of bands that 
can be observed in the spectra of the corresponding band 
systems. In those cases where there were not enough spec- 
troscopic data to determine woz or woyo, we did not deter- 
mine them theoretically. For the NO* molecule, we calcu- 


Table | 
\ 
tp. cm | | | &p 8p 
A 
1 108700 1 119000 33 
5 111700 1 420000 71 
3 114800 3 121000 54 
1 412100 5 122000 73 
3 113500 10 123000 103 
2 113000 16 424000 176 
7 113700 10 425000 293 
5 114100 3 426000 362 
8 415000 24 127000 45 
1 116700 428000 24 
3 417000 15 
5 117200 5 
3 418000 12 
A 
4 450000 45 490000 
2 155000 12 192000 
2 458000 20 || 193000 
14 160000 24 195000 
12 167000 2 4196000 
6 170000 14 499000 
6 173000 26 200000 
10 475000 6 205000 
48 480000 12 209000 
2 483000 10 213000 
184000 34 
12 186000 24 
18 187000 30 
= 0 5 99760 9 423000 
158,5 3 101150 15 424000 
45868 5 102800 29 127000 
43792 1 403870 9 429000 
33768 5 105300 48 129800 
76795 3 106700 57 131000 
86630 15 107500 48 432000 
88630 9 108200 134000 
“a 95480 5 108700 48 490000 
96930 3 113300 15 
97420 45 114000 37 
99090 15 116000 8 | 
or 
He 9 4 212700 2 254000 14 
26808 6 214200 6 255000 56 
26829 4 226800 2 4156000 138 
40467 6 228700 14 258000 14 
120000 12 231200 38 259000 14 
166000 10 233000 56 38 
185000 12 234000 12 266000 150 
189000 6 238700 12 268000 18 
195700 2 240400 6 276000 122 
203900 2 245900 26 292000 36 
207000 30 248200 16 
208400 42 251000 20 
212000 18 253000 28 
N 
4 96800 12 | 107400 4 
19223 6 97800 6 | 110000 14 
19231 4 99660 40 || 110200 30 
28840 6 403700 12 110300 48 
12 4104200 6 410500 26 
86130 2 4104600 6 | 112300 6 
86223 4 4104700 28 | 143000 102 
88110 6 104900 26 414000 100 
88153 4 105000 20 415000 92 
88173 2 405100 10 416000 340 
93582 2 406500 2 
94800 20 406800 20 
95500 12 407000 12 
eo 2 245700 6 334500 10 
474.5 4 267200 10 36300 12 
57192 2 287600 42 339800 44 
57252 4 297000 6 342000 10 
57333 6 301090 2 343000 32 
101030 10 310000 26 355000 32 
131000 2 311700 6 368600 12 
145880 2 314200 4 374000 16 
446000 4 317600 22 378000 56 
186800 4 18 387000 202 
203080 10 327000 2 
221300 2 330300 28 
230400 7 332800 22 
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at lated the quantities wx, for all four electronic states and energy of the electronic state is described by a Morse function 
nd we determined w, for the state *II. The constant w,2z, for in accordance with the formula 
oni a known dissociation energy of a given electron state was 
letermined by us under the assumption that the potential Were = (2/4D,) 11 
1. Molecular nitrogen N. Table 2 
DIN.) ==78715 (15 C),g-mole (chem}] (references 2 and 6) sid 
I(N.)==125702 cm~'=359283[gcal (15 C)/g-mole (chem)]. From I(N), D(N..) and I(N.*) 
Electronic state 
Name of constant 
| Bll, vz, N only 
0 49756 .5 59313.5 } 60000 68956 .6 63959 88984 
dt 2345.16 1446.46 1719.64 1515.5 1678.959 (1515.5) 2020 .00 
@oZo 14.445 13,929 14,47 11.5 13.3181 (11.5) 26.40 
—0.000509 0.000614 —0 002298 (—) 0.5542 
A 2.010 1.440 1.6380 1.48 1.637 (1.49) 1.8259 
a, 0.01869 0.013 0.0184 0.015 0.0224 (0.015) 0.0197 
dD, 5.8-10-* (5.61- 10-8) (—5.86-10-*) (5,56- 10-8) (6.14-10-*) (5.56-10-8) 6.0-10-6 
—0.001-10-§ (—0, 154-1078) (—0.0281 - 10-6) (—0.029- 108) (0.0645 - 108) (—0.029- 10-8) (—0.013- 10-8) 
‘im 50 26 30 30 25 30 4 
Source reference (2) (3) (2) (3) (3) (3) (3) 
2. Molecular oxygen O. Table 3 
D(O.)—41260 cm~'=117930[gcal (15 C) ‘g-mole (chem)] (reference 7) 
cm~*==282430[gcal (15 C) g-mole (chem)]. From I(O), D(O.) and D(O.*) 
a | Electronic state 
Name of constant — = 
0 7892 .39 13120. 908 36212.74 35780 .12 49363 
1568.33 1496.4 1418.7292 633.371 758.89 692 .08 
11.993 12.9 13.9662 17.2056 14.635 8.5653 
0.0517 —0,01075 —0.1205 --0.55 —0.37535 
1.445666 1.4264 1.400416 0.826 (0.902) 0.819 
ae 0.015791 0.0171 0.01817 0.02055 — (0.0245) 0.011 
D. 4.957-10-8 4.8652 5.356- 10-8 (10.731-10-8) -- (3.32-10-8) (4.38-10-) 
be —0.088. 10-8 —0.0104- 10-8 —0.077-10- (0.262- 10-8) (0.96- 10-8) —0.0784-10-¢ 
Vm 36 30 27 13 we 1! 21 
Source reference (12) (13) (12) (14) (14) (15) (7) 
3. Nitrous oxide NO Table 4 
D(NO)—52477 cm~'=149990 geal (15 C) g-mole (chem) (reference 7) 
I(NO)==76626.8 cm~'—219015 [gcal (15 C) g-mole (chem)] (reference 2) 
Electronic state 
60.6 441992 45440 52817 60862.8 
1889.88 2356 .82 1030.15 2304 2357.79 
13.972 14.48 7.458 23 15.87 
B, 1.7046- 1.9948 1.126 1.9912 1.9863 
| 0.0178 0.0184 0.0152 (0.0243) 0.0182 i 
| D, 5-10-8 5.6-10-8 2.65-10-8 5.8-10-* - 
be 0.02. 10-6 0.02. 10-8 (0.0405-10-*) (—0.063- 10-#) 
Ven ea 3 29 25 25 
| Source reference (8) (2) (3) (3) 
| 4. Molecular ionized nitrogen, N:” 
D(N.* )=70.358 cm~'—201098 gcal (15 C)/g-mole (chem) (reference 9) 
Electronic state 
Name of constant | (X#X+8) | Bert, | 
Ee 0 (50470) 9057 25566 (51618) 64547 
2191.02 (1023.5) 1887.93 2396.22 (1007.5) 2162.8 
13.19 (13.24) 14.91 24070 
B, —1_932 (1.315) 1.722 2.083 (1.433) 1.16 
ers 0.02 (0.0192) 0.018 0.0195 (0.0294) 0.05 
5.75 40-8 (7.427. 10-8) 4-40-* (6.15-10-*) (13.903- 10-6) (3.8- 10-8) 
0.29-10- (—0.0202-10-*) (—0,216-10-*) (0.216-10-*) (0,116-10-*) (0.6-10-*) 
<A 30 : 37 63 14 37 10 
Source reference (9) 9) (10) (9) (9) 


We estimated w, for *II from empirical relations. For certain 
electronic states, owing to interaction with the neighboring 
ones, it is in general impossible to select a series [7], with a 
small number of terms, suitable for all the vibration levels. 
In this case we select two series, one of which corresponds 
to arrangement of vibrational levels from v = 0 to a certain 
v,, and the other one from 0+: t0 tmx. This is equivalent to 
introducing a fictitious electronic state. Such fictitious 
states are set in parentheses in Tables 1 to7. The number of 
spectroscopically investigated vibrational levels for various 


tional levels above the investigated levels obeys the same 
law as that of the observed levels, then v,, can be determined 
from the relation 


= Do [12] 


Almost all the v,, listed in Tables 2 to 7 were obtained pre- 
cisely in this manner. However, Equation [12] had no real 
positive roots for certain electron states. In such cases 
Vm was determined from the following considerations. The 
dissociation energy of any electronic state should be equal 


electronic states is variable and fluctuates from one to several 
times 10. For the majority of the lowest electronic states 
usually we investigated a large number of bands, a fact favor- 
ing the reliability of calculations based on the spectroscopic 
data. However, the quantity vmx was not determined spec- 
troscopically for the majority of electron states of the mole- 
cules O2, Nz, NO, 02+, Net and NO*, which were of interest 
to us. If it is assumed that the arrangement of the vibra- 


to the sum of all the vibrational quanta AG(v) of the given 
kat! 
 D= D AG) 
v=0 


In other words, if we plot AG(v) vs. v, then the area bounded 
by the curve and the coordinate axes should be equal to the 
dissociation energy of the given electronic state. The solid 
curve in Fig. 1 corresponds to the investigated values of 
AG for the state X >; of the molecule N2+. If we continue 
this curve until it intersects the abscissa (dotted curve), then 
the area under the curve will be found to be less than 70,358 


cm~? (dissociation energy of the X state of the mole- 

: shyt % ah cule). This corresponds to the fact that Equation [12] does not 

have asatisfactorysolution. To obtain an area of 70,358 cm~' 

under the curve, we make a simple assumption concerning 

the course of the curve AG(v) above v = 30, namely, that 

BG) depends linearly on v, as indicated on the diagram by 

the dotted Jine. The constants wor and for the fictitious 
electronic state were determined from the equations 


G 


0 0 2 30 50 GO Wv 


- 1 Dependence of AG on the vibrational quantum number ‘ 
the state X* of the molecule N,* 


Name of quantity Numerical value 
5. Molecular ionized oxygen O: 1 light, c = 299,792 + 0.8 
D(O.,* }52283 cm~'==149.436 geal (15 C)/g-mole (chem) (reference 2) 2 Number of mole- 
Electronic state cules per g-mole 
Name of constant xn, > physical (1/g-mole) Non = (6.02472 0.00036) - 102% 
chemical = (6.02308 + 0.00033) - 
we wer Planck's constant, = (6.6252 0.0005)-10-" 
1859.87 1025.30 886 .6 1179.68 erg-sec 
18.58 1030 13.4 _17.08 con- k = (1.38042 + 0.00010)-10-"8 
2 § Electronmass,g = (9.1085 + 0.0006)-10-% 
2 (6.84-10-*) 4,88-40-* (5.33-40-#) 6 of the mass of 16g¢/(16 N.,) = 1.66028-10-* 
B, (—0.0356- 10-*) —0.095-10-* (—0.0955- 10-*) —0.0185-10-% j the oxygen atom, g 
constant, R = (8.31662 + 0.00038)-107 
Re, = (8.31436 + 0.00034) - 10? 
erg/mole-deg 
mimeo Gram-calorie,15C 1 g-cal (15 C) = 4.1855 abs joule 
9 Atmosphere, l atm = 1.01325-10* dyne/cm? 
physical 
Relations between energy unite 
0 37000 57000 74000 erg ev 
1 1322.4 
(1,96) (1.54) (1.35) (1.31) em~! 1.98620-10-% 1 
vr a, (0.0181) (0.174) (0.0196) (0.0163) he 
me (5.4-10-*) (5.3- 10-8) (5.5: (5.5- 10-4) = 1.43883 em-deg 
(73) (39) (42) @) 
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AG = height of the first vibrational saa above the zero level 
(HF on the diagram) 
Dy = area of the triangle FBH 


Not more than four rotational constants are usually deter- 
mined spectroscopically: B., a, D, and 6,. Where there 
were no data on these constants, they were calculated from 


the 


_ 3 [14] 
« 
Swetre 


In most cases the formulas give the correct order of the 
constants. The rotational constants B, and D, cannot always 
be well represented by Equation [6]. Since the deviations 
from the ordinary dependence of B, and D, on » correspond 
to the same deviations for the vibrational constants, the 
introduction of the fictitious electronic states results simul- 
taneously in convenient methods of calculating both the rota- 
tional and vibrational parts of the statistical sums. 

As already mentioned, the rotational parts of the statistical 
sums of the molecules were calculated in this investigation for 
a finite number of values of J. The maximum value of J at 
which the molecule is still stable was estimated in the follow- 
ing manner. 

The potential-energy function of the rotating molecule 
can be represented by 


h 
U(r) = DA — + +1) [16a] 


= (r—1)/te 
= equilibrium state between the nuclei 
and 


We Wels 


U(r) has a maximum located above the dissociation limit. 
This maximum forms a potential barrier, beyond which there 
exist stable molecular states with energies greater than the 
dissociation energy. The position of the maximum for each 
value of J(J + 1) can be found from the condition 


which after transformation reduces to 


From the value of rmx thus obtained for the given J(J + 1) 
we can calculate 


Um = — e~28rm—re/re) 4. +1) [17a] 
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h 
where al 


1\2 


> 
B.J(J +1) +1)? +... [18] 


we can obtain the value of J(J + 1) that corresponds to the 
given v. Equations [17 and 18] were solved graphically. 

The scantiest data on spectroscopic constants were those 
for NOt. 

Not one of the electronic states of NO* has yet been investi- 
gated spectroscopically. There are only data on the location 
of the first four electronic states, obtained from a spectro- 
scopic investigation of the Rydberg series of NO*. For 
three out of these four electronic states, the existence of which 
we have verified from the convergence of the Rydberg series 
of NO+, we know the values of w., also obtained from spectra 
of NO. For the lower electronic state we know also the 
equilibrium internuclear distance 


= 1.073 A [18a] 


The method of electronic impact was used to determine the 
ionization energy of NO, found to be 9.5 ev. From this, 
using the equation 


I(NO) + D(NO*) = D(NO) + I(O) [18b ] 


we can obtain 


D(NO*) = D(NO) + I1(0) — I(NO) = 10.5ev_ [18c] 


All this disjointed information was gathered together by 
Rosen (3). The calculation of the statistical sums of NOt 
on the basis of such information may lead to an error, due to 
lack of knowledge of the electronic terms, by factors of two, — 
three, or four. We therefore attempted first to determine 
the types of electronic terms of the known four electronic 
states. One method of determining the possible types of 
electronic terms is to construct the electronic configurations 
by successive addition of the electrons. Such a compilation 
of the electronic configurations for the molecules can be 
carried out in two limiting cases; first, when the distance 
between nuclei is small compared with the distance between 
the electron and the nucleus, and second when the distance 
between nuclei is large compared with the distances between 

- the electrons and the nucleus. 

For the atoms N and O+ “joined” together (r, ~ 0), we 
would obtain the following sequence of filling of the electron 
shells 


[18d] 


For the N and O+ atoms far apart we would obtain the 
following sequence of filling of the electron shells 


6280 )?( 02 {18e] 


The order of arrangement of the electrons in an actual 
NO?+ molecule will be intermediate between the system in 
which the internuclear distance approaches zero, and in a 
system in which it goes to infinity. The NO* system is 
isoelectronic relative to the molecules Nz and CO, which have 
14electronseach. For the lower electronic state of co(?>>*) 
the value of r, is 1.1289 A, and for N,(? , Pe it equals 1.094 
A. For the lower electronic state of NO*, as already men- 
tioned, r, = 1.073 A. This gives us the right to assume that 
the lower electronic configuration of NO+ will coincide with 
the lower electronic configurations of CO and N2 


An important argument in favor of such a lower electronic 
configuration of NO* is the fact that it corresponds to the 
removal of one upper electron m2py from the NO molecule, 
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which has a lower electronic configuration ges “ayeethiute 
o2pn )*(42pn) [18g] 


This method of determining the electronic configuration 
(we have in mind the method of removing the electron from 
the electronic ground state of the un-ionized molecule) is 
found to be valid in examples of other molecules (for example, 
determination of the electronic states of N2*). 

The first excited electronic configuration of NO* may be, 


apparently 
) {18h 


or 


KK(02s0)*( 02s 72po)( o2pn) [18i] 


Removal of one of the o2po electrons of NO results in the 
electronic configuration 


(18j] 


which coincides with the first of the possible excited electronic 


configurations of NO*+. One can therefore assume the 
electronic configuration 
[18k] 


to be the first excited electronic configuration of NO*. 


es The next excited electronic configuration of NO*+ may be 


I 


or 


Removal of the (#2po) electron of the NO molecule leads to 
an electron configuration 


KK(o2s0)*( 


which coincides with the second of the possible ae “a 
tronic configurations of NO* from the series (I, II, III). 
We thus conclude that the possible electron configurations 
of NO* are: 


Normal state 


Second excitation 


Corresponding to these electronic configurations are the 
following terms: 


Normal state 
First excitation #II; *II 


The possible electronic terms can also be obtained from the 
products of the dissociation of the NO* molecule, by using 
the Wigner-Witmer moment-summation rules. In Table 9 
we list the possible products of dissociation of NO+, the 
possible electronic terms corresponding to them, and the 
relative arrangement of the limits of dissociation. 

Comparing Table 9 with the possible terms obtained from 
the electronic configurations, we can state that the term of 
the normal state of NO* is the term *}>+. Of the remaining 


terms of_the first limit of dissociation 


te 


Height of the 
products “the limit of the lower, 
electronic state, cm 
N(#S) + O* (4S) 1; 3; 55 0 
N (2P) + 0+ (4S) 3; 55+. 35 19223 
N+P) +- (2); “Et; “I (3); “A (2); 23376 


the only stable one is apparently > Par since it has only one 

pair of uncompensated spins. The terms *)>+ and 7)>+, as 

in the case of the Nz molecule, are unstable, since the number 
of uncompensated spins predominates in them over the 
number of compensated spins. We are led to the same result 
also by the fact that the terms *}>+ and ’}°>*+ are not included 
among the terms obtained on the basis of comparison of the 
electronic configurations. The term *)°>+ is encountered 
in Table 9 in all four cases. However, to satisfy the rule of 
noncrossing of the electronic terms, it is necessary to assume 
that term *}°+ of the second excited electronic state has 
dissociation products 


N(4S) + O+(48) 


The term of the first excited electronic state *II is en- 
countered in Table 9 for the dissociation products oy 


N*+(*P) + O(*P) 
N(?2D) + O*(48) 
N+(®P) + O('E) 

Using the noncrossing rule, it can be said that ‘II has dis- 
sociation products N+(*P) + O(’P). The same dissociation 
products are possessed, apparently, also by the term 'II, 
since it is not encountered in any of the rows of Table 9 
(with the exception of the second one). 

We can thus apparently state the following. Located far 
above the ground state 1}°+ of the NO* molecule are the 
terms "II, *II and *}>+, which are relatively close to each 
other. This is connected with the fact that the orbital 


o - No 
o- 
x - NO 
o - NO* 
U - 
0 


[A] 


Fig.2 Dependence of r, on we for Nz, N.+, NO, +, and O,+ 
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(x*Px) is quite far above the closed electronic configuration 
of NO+. The orbitals o2Po and r2Po, from which the 
terms ‘II, *II and *}°* are obtained, lie close to each other. 
Using the Hund rule, by which the state with the greatest 
multiplicity lies below the remaining ones, we can say that 
II lies below 4II. The state is apparently located 
above *II. It cannot be located below ‘II, since it is obtained 
by transfer of a deeper electron from the orbital r2Po to the 
orbital 72Py, whereas *II is obtained by transfer of an elec- 
tron from the orbital o?Po to the orbital 72Px. This is 
supported also by the fact that r, for NO is less than r, for 
CO for the lower electronic state; i.e., the molecule NO is 
closer than CO to the limiting case r, ~ 0, and with this the 
orbital ¢2Po becomes more and more remote from 72Po. 
The term ‘II is apparently located above *II. The only proof 
in favor of such a relative placement of *}°* and ‘m1 is the 
analogous placement of the states *}°+ and 'II in the iso- 
lectronic molecules CO and Ne. It must also be noted that 
the states *}°+ and 'II are high (respectively 57,000 and 
71,000 em-!), and an error in the choice of the term results 
in either a factor 3 or a factor 2 in front of the statistical 
sum of the electronic state with 7 = 57,000 em-. 

To calculate the rotational] part of the statistical sum it is 
necessary to have at least values of r. from which to calculate 
the quantity B,. For NO+, we know r, for the state '}°*. 
An estimate of the values of r, for the three other states was 
made from empirical laws. For the molecules No, Not, Oo, 
O.* and NO we noted the following regularity. The de- 
pendence of r, on w, for the electronic states located not 
more than 70,000 cm above the ground state, is a relatively 
smooth descending curve. With this, all the values of r, fall 
on one curve for Nz and Net, and on another parallel curve 
for O. and O,+. The curve for NO lies between the No, 
N2* and Os, O.* curves. One known value of r, for NO fits 
satisfactorily on the curve for NO. From the curve for NO 
we determined all the missing values of r, for NO+. Fig. 2 
shows all the curves and points from which the values of 7, 
for NO* were taken. 


3 Effect of Ionization, Coulomb Interaction and 
Degeneracy of Electron Gas on Values of 
Thermodynamic Functions 


We shall show below that it is possible to neglect the 
Coulomb interaction between the charged particles when 
estimating the degree of ionization. Let us denote by No, 
Ni, No, ..., Na,...,Nm the number of neutral, singly-, 
doubly-, ..., n-fold ionized atoms per gram-mole of the 
gas; we denote by N the total number of atoms in a gram- 


mole of gas = 


by N, the number of electrons in a gram-mole of gas, and by 


The ionization process can be considered as a system of 
chemical reactions 
(1) 
At+ At++ (2) 
n( +) (n—1)( +) (n—1)(+) 
A +I[A ] (n) 
AmMH 4 4 (m) 


The condition of thermodynamic equilibrium for such a sys- 
tem of reactions is expressed in the form of a system of alge- 
braic equations 


p 
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considerably simpler. 


where Co, Ci,...,C,,..., Cm are concentrations of the 
neutral, singly-, doubly-, and m-fold ionized atoms in a gram- 
mole of gas and 


C 
K p™ 


concentration of the electrons 

equilibrium constants for the corresponding reac- 
tion, depending on the potentials of the ionization 
and the temperature 

Pp = pressure 


Q 
| 


NN +N) 
+ N,) 
N,/(N + N,) 


The m equations of equilibrium and the m + 1 unknown 
must be supplemented by the condition 


2 


On +C 


and the law of conservation of charge a a 


which expressed in terms of the concentrations becomes 


Transforming from concentrations to degrees of ionization, 
Q, Q1,..., Q@m, we obtain the following system of m + 1 
equations for m + 1 unknowns, 


+ 2a%2.+.. 
(1 + a, + 2a2+... 


+ nQn+.. 


+ Mam) _ 


MAm)An-1 


am 
1 


A calculation of the equilibrium constants shows that 


K, (n) > (n+1) 


and the ionization process proceeds, therefore, in steps as the 
temperature is increased. First the gas is practically all 
singly ionized, then it is all doubly ionized, and soon. With 
such an ionization process, direct calculations show that the 
system of m + 1 equations for ao, a1, ..., Qn, ... &m becomes 
When the (n — 1)th degree of ioniza- 
tion is complete and the nth degree of ionization begins, the 
following relations hold 


3. 

An + Qn-1 = 1 

Onti = Anpo = =... = Am = O 


With this, the (m + 1)th equation for a becomes a apenne wat 


two equations for a,_: and &, 


‘ee, 
An[(n — + _ Ky 

+ On = 1 


an 


n-1 n—1\? nKk,™ 


for 

Dove 

wer 

> 

ac 
as 

ber 

the 

ult 

led 

| 

the 

has | 

a 
| 
J 


When n = 1, the expression obtained goes into the well- 
known Saha formula 


1 


a = 
Vi + 


7 Table 10 gives the results of the estimates for a, a, and 
a. The estimates were made for nitrogen. The ionization 
potentials of N, O and A as well as the ionization potentials 
of N+, O*+ and At are approximately equal. Therefore 
Table 10 is suitable also for estimating the ionization of O 
and A. 

The deviation of the gas from an ideal gas, due to Coulomb 
interaction, is characterized in the Debye-Huckel approxi- 


ionized gas 


kT 
= positive or negative numbers, characterizing the charge 
of the corresponding ion 
number of particles of the nth sort in a gram-mole of 
gas; the index n is not connected with the multiplicity 
of ionization 


where 


2n 


The results of the calculations of the corrections are listed 
in Table 11. 

It is seen from Table 11 that up to temperatures of approxi- 
mately 12,000 to 14,000 K one can neglect the Coulomb inter- 
action in the pressure interval from 0.001 to 1000 atm. The 
degree of degeneracy of the electron gas is characterized by 
the parameter 


Table 10 
0.001 0.01 01 1 10 100 300 1000 
0.16 0 
1.4] 0 
97 0} 78| 0} 36] 0 |12 4 
100 0} 97| 78] 0 12 
99 1/100] 95] 0 
86 | 98] 0 
45 | 55| 88] 12] 99] |95.6|0 65/35/86 
13 | 87] 57] 43] 92] 8|99 | 4/100 |o 
3 | 97] 16] 84 | 67| 33 |96 | 4 99.50.5 38 | 62 
Ap 
A 
12 000 0.001 0.7-10-* 0.07% 
14 000 0.01 13.6-10-6 0.14% 
Table 12 
| | | | | 
0 Ar 0 N;* 359 
0 Nt 447 889 O,* 282 430 
1 467 o+ 372 901 NO+ 219015 
112 492 Art 363 307 
58 965 e 


mation by a correction term in the equation of state eal the | 


At T = 25,000 K and p = 1000 atm we have D = 0.003< 
1; ie., we can neglect the degeneracy of the electron gas. 


4 Formulas for Calculation of Enthalpy, 
Specific Heat and Equilibrium Constants in 
Temperature Interval from 1000 to 12,000 K 

It was shown in Section 3 that at these temperatures the 


gas can be considered ideal, and therefore the formulas for 
the calculation of the enthalpy and heat content have the 


form Tt’ Sa 
2 
> 


where R is the gas constant, and the values of the constant :C 
are given in Table 12, in g-cal (15 C)/g-mole (chem). 
The formulas for the equilibrium constants are of the form: 


For the reaction 2N = + D(N2) 
Qn? 


Ne 


the reaction 20 = + D(O2) 


K, = 0.474542T*/t ¢-113,257/T atm 


= 0.5 79281 7°/2 Qo? 59,366/7 atm 
For the reaction N + O = NO + D(NO) we a 
Qno 
For the reaction of ionization A+ + +1(A) 


K, = 6.57936-10~7 atm 
A 


where 
I(N) = 168,840 deg 
I(O) = 157,036 deg 
I(A) = 182,890 deg 
I(Ne) = 180,864 deg 
I(O2) = 142,176 deg 
I(NO) = 110,253deg 
pal 
eu) 
5 Solution of Equations of Chemical © wif 
Equilibrium of Ideal Gas Systems 


For a thermodynamic analysis of reacting gas systems it is 
necessary to calculate first the composition of these systems. 
The principal solution of this problem is known. However, 
in the case of a large number of reactions and components, 
a direct solution of the corresponding equations leads to cum- 
bersome derivations which are difficult to visualize. 

In this section of the article we discuss briefly a method of 
- symmetrizing the equations of chemical equilibrium with 
and without allowance for ionization 

It is known that the closed system of equations of chemical 
equilibrium (ECE) of an ideal gas system consisting of non- 
ionized? components, containing a different types of atoms, 
consists of one equation of Dalton’s law (EDL), a — 1 equa- 
tions of material balance (EMB), and r = k — a equations 
using the law of mass action (LMA). Here r is the number 
_of reactions that proceed in this system. 


ae 2 The case of ionized components will be considered later. 


ARS JourNnaL SuPPLEMENT 


- 
‘ 
> 
— 
T° 
1 
1 
2 
2 
2 
> 
a 
> 
3 
» 
1 
* 


< 


he 
or 
he 


We introduce the symbols: 


K; = component (i = 
An = atom m (m = l,. 
Nim = number of A» in the molecule K; 
v,; = stoichiometric coefficient of K; in the reaction 
a 
The molecule K; has the form >> nimAm, while the reac- 
mal 
k phe 
tions — TheECEhavetheform 
i=l 
2=1 (EDL) ad [19a] 
i=l 
mints | Yon Niti = Pmi (EMB) [19b] 
xi ms (LMA) ‘[19¢] 


pu = 1 


The physically obvious uniqueness of the solution of the 
ECE was proved in (15). 

Equation [19] is difficult to solve because of the non- 
linearity of Equation [19c] and the asymmetry of Equations 
[19a and 19b]. 

The principal numerical methods used to solve the ECE are 
the iteration method and the Newton-Rafson method. The 
need for using numerical methods in solving Equation [19] 
arises from the nonlinearity of Equation [19c]. We iterate 
here the a components whose content predominates in the 
region of the p-7' diagram where the solution of Equation 
[19] is sought (we have in mind here the values of z;). Differ- 
ent regions of the p-T diagram have different sets of such com- 
ponents. 

We shall henceforth segregate the set of the so-called pre- 
dominating components, the number of which equals the 
number of various types of atoms, i.e., a. This set is as- 
sumed ordered in such a way that the mth component of the 
set has the maximum content of the mth type of atom, 
compared with all other components. Under this definition, 
this set may contain components the absolute content of 
which is small. We shall assume that the predominating 
components are the first a components from among a total 
of k components. The asymmetry of [19a and 19b] is re- 
moved by solving them with respect to the predominating 
components, i.e., by obtaining symmetrical equations for the 
concentrations of the predominating component (EPC). 
For this purpose we introduce the concentration of the pre- 
dominating components 2,,° in the absence of the remaining 
components, where 2z,,° are determined from [19a] or [19b]: 


0 
= = pm (pu = 1) 
From this, obviously 


We furthermore introduce 
aye 


t=1 


In view of the foregoing definition of the predominating 
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components, Nmm * 0. Then [19a and 19b] assume the form 


m=1 


1+, 
From this we readily obtain the EPC that are equivalent 


to the EDL and the EMI a - 


a k a 


The EPC becomes simple if each K,, consists of only Am. 
Then nin = 0 (l ¥ m;1,m = 1,..., a) and Equations [21b 
and 22b] become 


mm > [21b’ 


m=1 l=a+l1 m=1 mm 


In this case the EPC express the concentrations of the 
predominating components in terms of the concentrations of 
the remaining components. Since in most applications the 
LMA make it possible to express the concentrations of the 
remaining components in terms of the concentrations of the 
predominating components, the systems [22a and 19c] are 
suitable for generalized iteration. We note that it is always” 
possible, by combining Equations [19c], to recast the LMA 
in the most symmetrical form, which permits a simple ex- 
pression of the remaining components in the term of the pre- 
dominating ones. 

The calculation of the iterations, which is known to be 
purely algebraic, is simple, but since the iterations converge 
only in the centers of the regions where the components 
predominate on the p-7' diagram, their application is limited 
in character. The Newton-Rafson method is much more 
promising. 

The principal advantage of this method is that the region 
of its applicability goes far beyond the boundaries of that 
region of the p-7’ diagram where the components, relative 
to which the EPC are solved, predominate. This is con- 
nected with the differential character of the method. The 
amount of calculation necessitated by this method is con- 
siderably less than that involved in simple calculations based 
on the EDL, ENI and the unreduced LMA, thanks to the 
symmetry of the ETC and of the reduced LMA, provided 
the components are suitably numbered and the computation 
program is suitably detailed. Finally, the Newton-Rafson 
method gives information necessary for the calculation of 
the specific heats, the velocity of sound and other thermo- 
dynamic functions of the reacting gas systems in the equi- 
librium state. 

We shall detail in the following a program for the calculation 
of the ionizations of arbitrary orders of all & components. 
Bearing in mind that the composition is independent of the 
specific choice of the complete system of independent reac- 
tions, we can write the equations of the ionization equilibrium, 
with allowance for n-fold ionization, as vis; 


~ 
f= Pp M1 


1 


| 
Ss. 
r, 
8, 
Ss, 
er 


where z, is the concentration of the electrons 
n k 
t= i=1 
Equation [23b] expresses the law of conservation of electron 


charge. 
Multiplying (23b] by z- and making use of [23a], we ob- 


[23b ] 


Bi = molecular weight of the component 7 Ges 

N; = number of moles of the component 7 

Xi = N,/N, the concentration of the component 7 
H;,U; = molar enthalpy and internal energy of the same 


component 


The symmetrical equations of chemical equilibrium (ECE) 
consist of symmetrical equations of the predominating com- 
ponents (EPC) and the reduced equations of the laws of 


(5), 


_ number of moles in the system), we get 


tain 
k mass action (LMA); see Section 5 
where Kp, = I = 1,...,R) (EEML) [27b| 
If a:,..., are known, then for specified values of Kp, te 
Equation [24a] is an algebraic equation of degree n — 1 a = numberof different atomsofthesystem a. 
in z.. It is generally solved by numerical methods. For R = number of independent reactions, of which the first r are 


ordinary reactions, and the remaining ones are 


ionization reactions 


the first-order ionization we have 


After determining z,., the ionized components are found from 


24¢ 
_ If we go over from z; to N; and use the equation of state 


pV = NRT, then Equations [27] assume the following form 


the _ Nn + Vn = constant [28: 


The program ra ions the a ECR is as follows: Here W,, stands for ®,, when z; is replaced by N 


We specify the predominating components 2, ... , Ya; the From rae [25-28] we obtain 
transformed Equations [19c] are used to find ra+1,... , 
Equation [24a] is used to find and Equation [23a] is used Coit: + + Hy; + HE 
to find r4+1,..., We thenapply the Newton-Rafson i=1 
method to Equation [22a], obtaining new values of ,..., Za, h 
and so on until the necessary accuracy is attained. — he we a 
kK 
ame 
6 Thermodynamic Analysis of Reacting Ideal 2» N 
Gas Systems in Equilibrium State lt 
Let us consider first the calculation of molecular and specific 
heats C, and C, for systems of reacting ideal gases in equilib- gure 
rium. 
The calculation of these quantities is discussed in (16). It C= xX Cots — R + 2 Hw; — RTZ a . 
be shown below that this calculation is exact. inti 
The formulas for the specific heats c, and c, have the form y 1 (2 7) - > 


Thus, the calculation of C, and C, reduces to ve y; and 
v;, respectively. Taking logarithms and then differentiating 


(25a ] 


where h and wu are the specific enthalpy and internal 
energy of the system. Multiplying and dividing Equation 
[25a] by the constant mass of the system M = uN (u 


> uit; = molecular weight of the system, N = > N; = 


i=1 i=] 


Equation [27b] with respect to 7, with p held constant, we 
obtain 


The bars denote the substitution of y; for z;. In Equation 
K K nie a [29b] we use the well-known relation 
where 
K = number of components, of which the first k are the where Qe. = »» viH; is the heat of the reaction s. 
un-ionized components Analogously, we obtain from [28] at V = constant 
’ The molecular quantities are denoted by capital letters and 
tm + Vn =0 [30a } 


the specific by corresponding lower case letters. 
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: 4 
- 


K 
MU ape = Ly = Qays — [30c] 
i=l 
Here the bars denote the substitution of v; for z;. A direct 


olution of systems [29 and 30] is not advisable. It is simpler 
) first find y, from Equation [29a] and v, from [38], and find 
from these the values of yay, ..., yx from [29b] and of 
41, . . +», Ux from Equation [30b]. Such a sequence of 
eration (from the EPC to the LMA), which is natural 
for the forms [29 and 30], is rational whena <R. Whena> 
R we should use the reverse sequence of operation (from the 
EECM to the EPC). Actually, at a > R one can confine 
one’s self to a solution of only the LMA, transformed in 
suitable manner. The case a > R will be considered some- 
hat later. 
Retur ing to Equations [29 and 30], we note the sub- 
antial obvious circumstance: The matrix of the solution 
ised on Equation [29a] coincides with the matrix of the last 
ep in the determination of dz, by the Newton-Rafson 
ethod by means of Equation [27a]. 
To calculate sas when a > R we begin with the system 


7 


Na t+ V,, constant [31a] 
i=1 


which differs from Equation [30] by the obvious substitution 
and the replacement of the index of K,. We note that the 
conditions of conservation of the number of atoms of each 
kind must be observed for each reaction (see Section 5) 


K 
>, Mittin = 0 [32] 


To calculate C, from Equation [26] it is necessary to deter- 
mine u,; = (1/N) (ON;,/OT),. After taking the logarithm of 
31b], differentiation of Equation [31] with respect to 7’, with 
p held constant, yields a system of linear equations for u; 


K 
= 0 


é (33a ] 
If we furthermore represent u; in the form 


then Equations [33] are automatically satisfied by virtue of 
Equation [82]. 
To determine IT’,, we derive the system of equations 


here 


[36] 


To calculate C, for a > R we start, as before, with Equa- 
tions [28], that is to say with Equations [30], which are written 


in transformed form, in analogy with Equations [33] 


Nims = 0 


[87a] 
AU, v; 

RT? ee > 


Here, too, the index of AU has been changed. 
Further operations are analogous to the calculation of C,, 
and we obtain for the representation of 


equations that permit us to determine 
R 
=) AU, 
where 
1 
—~=),“*=— +4», 39 
Set p> Fs: 
From Equation [25] we obtain for C, 
= >, Cun + AUN [39b] 


i=1 


Epstein (16), who proposed a system analogous to the 
foregoing for a > R for the calculation of C, and C,, did not 
obtain the systems [35 and 39], inasmuch as he implied ma- 
terial independence, which in general is incorrect. 

The term “materially independent reactions” denotes 
reactions for which the corresponding sets of components 
do not overlap. (Each reaction has its own corresponding 
set of components, which participate in this reaction.) 

Epstein determined I’, and , from the equations _ 


Q, 
[39¢] 
without solving, thus, the systems [35 and 39]. If the inter- 
action of the reactions is weak, Epstein’s calculations are 
sufficiently accurate. This occurs when the regions of the 
maximum intensities of all the reactions do not overlap on 
the p-T diagram. For a single reaction in this system, nat- 
urally, Epstein’s calculation is correct. It is therefore no 
wonder that theory and experiment are in good agreement 
in the experiments of McCollum (17), where only the reac- 
tion 2NO.-N.O0 = O was investigated. As to the strongly 
interacting reactions, naturally, Epstein’s calculations are 
not valid for them. In the calculation of the specific system 
of reactions discussed in this paper, we found a discrepancy 
up to 40 per cent between the values of the specific heats 
calculated by the proposed scheme and those calculated by 
Epstein’s scheme, over a sufficiently large region of the p-T 
diagram. 

We note that it follows from Equation [35a] that the ma- 
terial independence of the reaction is not a sufficient condi- 
tion for the correctness of Epstein’s scheme in the calcu- 
lation of C,. In fact, this calculation leads only to v,.0,; = 0 
(s t). It is obviously necessary to have also the condition 
that the number of particles be constant during the occurrence 
of all the reactions, with the exception of one. This last 
condition reduces to v, = 0, where s runs through all the 
values from 1 to R, with the exception of a single one. The 
material independence of the reactions is sufficient for 
Epstein’s scheme to be correct only in the calculation of 
C., as can be seen from Equation [39a]. 
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From Equation [25] we obtain for C,, the following formula 
K R 
30a 


T 
+ N (= ), 


RT 


The value of the velocity of sound, calculated in accordance 
with the proposed scheme, differs from the value given by 
Epstein’s scheme, although this difference is less than that 
obtained in the calculation of the specific heats. 

Let us turn now to an examination of other thermodynamic 
characteristics. 

It is easy to see that the solutions of the systems [29 and 
30] make it possible to find partial derivatives with respect 
to T, with ¢ or V held constant, for all quantities that char- 
acterize the investigated system. We shall confine ourselves 
to a mention of the following equations pel 


oh 
(2) 


oh R 
= 
+ Ps (1 + T2,) 
op p 
2) T (1 + ent 


Analogously, we can obtain (0s/07’)»,». 

If derivatives with respect to T of H; of orders higher than 
the first are specified, it is possible to find for the quantities 
partial derivatives with respect to T (with p or V constant) 
of higher orders, but not exceeding the order of the derivative 
of H; with respect to 7. Corresponding systems can be ob- 
tained, by differentiation of suitable order, of the systems 
[29 and 30]. 

We can find equally well the partial derivatives of arbi- 
trary quantities with respect to p or p for 7 = constant. 
For example, to investigate the Joule-Thomson effect, it is 
essential to know (O0u/Op)z7, the equation for whichis 


i=1 


ll 
Rie 


The corresponding system for v; differs little in structure 
from the system [30], namely: All the Q.+, vanish identically 
and T is replaced by p; i.e., it is possible to calculate (0u/Op)r 
by using the ready-made program for solving the system [30]. 

Finally, using the properties of Jacobians (18), we can cal- 
culate the first-order partial derivatives of any arbitrary 
quantity (f) with respect to another arbitrary quantity (zx) 
for constant third derivatives of quantity (y). 

We note that the partial derivatives of any quantities with 
respect to arbitrary quantities at T = constant can be found 
for an order as high as desired. 

On the basis of the procedure expounded in the previous 
sections for the calculation of the thermodynamic properties 
of a mixture of gases capable of reacting chemically and ca- 
pable of being ionized, the authors have calculated the thermo- 
dynamic properties and the composition of atmospheric air 
in the temperature range from 1000 to 12,000 K and in the 
pressure range from 0.001 to 1000 atm. 

The volumetric composition of air under normal conditions 
was assumed to be as follows: Nz—78.08 per cent, O.—20.95 
per cent, with all the other components replaced by argon. 


7 Discussion of Results 


The formula for the velocity of sound is(16) 


u—internal energy (g-cal/g) 


The error in such an approximation is offset by the variation 
of the composition of air with altitude. 

Preliminary estimates and subsequent exact calculations 
have shown the following: 

1 Up to 6000 K, at pressures from 0.001 to 1000 atm, we 
can neglect the effect of ionization on the thermodynamic 
properties and the composition of the air. 

2 From 6000 to 12,000 K it is possible to neglect the 
effect of ionization other than that of first order. In connec- 
tion with this, the temperature interval from 1000 to 12,000 K 
is broken up into two stages: From 1000 to 6000 K withi- 
out allowance of ionization; from 6000 to 12,000 K with al- 
lowance for single ionization of all the components. 

The calculation relative to T is carried out every 100 deg Kx 
(a total of 111 values of 7’). The calculation with respect to 
p is in accordance with the following scale: 


32 rs From 1 X 10k to 2 X 10k—every 0.2 X 10k hindi 


From 2 X 10k to 5 X 10k—every 0.5 X 10k sats 
From 5 X 10k to 10 X 10k—every1 X 10K 
where k = —3, —2, —1, 0, 1 and 2 (a total of 16 points fc 
each order of p—a total of 97 values of p for each temper 
ature). A total of 10,767 points was calculated. 

In the temperature interval from 1000 to 6000 K six values 
of the molar fractions of Nz, O2, NO, N, O and A were calc 
lated at each point. In the interval of temperatures fror 
6000 to 12,000 K were calculated 13 values of molar fractior 
of Nz, O2., NO, N, O, A, O.t, NOt, N+, At ande. I 
the interval of temperatures from 1000 to 12,000 K, the fo 
lowing functions were calculated at each point, in ‘additio 
tothecompositionoftheair: = 


u—molecular weight (g/g-mole) 
c,>—specific heat for constant pressure (g-ca]l/g-deg) 
Cy—specific heat at constant volume (g-cal/g-deg) 

elocity of sound (m/sec) 
p—density (g/cm) 


The results of the calculations were used to plot 44 curves 
of which we include in this article curves for xn,, Yo, 2n¢ 
Zn, Lo, La, Le, h, Cp, a, p and s (Figs. 3-10) as functions of the 
temperature for pressures 0.001, 1 and 1000 atm. 

An analysis of the composition of the air based on thes 
curves shows that an increase in the temperature increases 
the content of first the atomic and then the ionized com- 
ponents; an increase in pressure, on the contrary, leads to a 
reduction in the contents of the ionized and atomic compo 
nents and to an increase of the molecular components. Fo 
all pressures, an increase in temperature develops first th« 
reaction of dissociation of molecular oxygen. The dissoci 
ation of Nz begins considerably later, since N2 has almost 
twice the dissociation energy of O2. The reaction of produc- 
tion of NO occurs simultaneously with the reaction of dissoci- 
ation of O2. Usually the content of NO in air reaches, with 
increasing temperature, a maximum value during the time 
when the dissociation reaction of O2 is the most vigorous, after 
which the content of NO decreases with increasing temper- 
ature. 

The ionization of atomic and molecular components begins 
in practice with 6000 K. As the temperature is increased, 
the percentage of ionized components increases, and at p = 
0.001 atm and T = 12,000 K the air is found to be fully 
ionized. 

The enthalpy and entropy curves increase monotonically 
with 7 and decrease monotonically with increasing p. It is 
possible to segregate on the enthalpy curves sections where 
the enthalpy increases markedly. The first such increase, 
measuring from the beginning of the curve, is due to the dis- 
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Fig. 3 Molar fractions of molecular and atomic components of 

ir and electrons as a function of temperature at p = 0.001 atm 
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Fig. 4 Molar fractions of molecular and atomic components of 


air and electrons as a function of the temperature at p = 1 atm 


1000 4000 7000 10000 12000 

Fig. 5 Molar fractions of molecular and atomic components of 

air and electrons as a function of the temperature at p = 1000 atm 
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Fig.6 Enthalpy of air in gcal/g as a function of the temperature 
at p = 0.001, 1 and 1000 atm. For comparison, the value of h 
calculated from Equation [41] is shown 
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sociation of N:, and the second is due to the ionization of the 
components of the air. 

The curves for the specific heats have several maxima. 
The first maximum from the beginning of the curve corre- 
sponds to the dissociation of O2, the second to that of Nz and 
the third to ionization. 

The molecular weight of the air diminishes with increasing 
temperature. The curves for u have sections where they de- 
crease rapidly. The first section from the start of the curve 
is due to dissociation of O2, the second to that of N2 and the 
third to ionization of the components of the air. 

The curve for the density is not reproduced here, in view 
of its practical] linear dependence on the pressure. 

The curves for enthalpy, specific heat, velocity of sound 
and molecular weight show also for comparison the values 
of h, cp, a and yp calculated in accordance with the following 
formulas 


RT 
a= 1.4 = 28.966 [41] 


These values of the thermodynamic functions, as follow for 
example from the curves, differ considerably from the true 
values. 

The authors have calculated the values of cp, ce, y and a for 
25 points using the Epstein scheme. Table 13 lists the re- 
sults of the comparison of the calculations based on Epstein’s 
schemes and on the scheme used by the authors. The com- 
parison is given for c, and c,. The differences in y and a, 
calculated in accordance with the two schemes, can be 
neglected in practice, something that cannot be done obviously 
as regards toc, and cy. 

The mathematical accuracy of the calculations of the 
thermodynamic functions and of the composition of the air 
is quite high and exceeds the physical accuracy. The sources 
of physical errors are inaccuracies in the calculation of the 
thermodynamic functions of the components of the air, the 
approximate nature of the determination of the composition 
of the air, mentioned at the beginning of this section, and 
also a certain incorrectness in the Clapeyron equation at low 
temperatures and high pressures (owing to corrections of van 
der Waals kind) and at high temperatures (owing to correc- 
tions of the Debye-Hiickel kind, which take into account 
the Coulomb interaction of the ionized particles and elec- 
trons). 

The computations were performed on a high speed elec- 
tronic computer, constructed at the Institute of Precision 
Mechanics and Computer Engineering of the USSR Academy 
of Sciences. 
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Fig. 7 Specific heat of air at constant pressure in gcal/g-deg 

as a function of the temperature at p = 0.001, 1 and 1000 atm. 

The value of cp calculated from Equation [41] is shown for 
comparison 
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Fig. 8 Velocity of sound in m per sec as a function of the 

temperature at p = 0.001, 1 and 1000 atm. The values of the 

velocity of sound calculated from Equation [41] are shown for 
comparison 
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10.1 
oe The results of the calculations were published, in the form 
— of tables of thermodynamic functions and of the compositior 
Petia of air in the interval of temperatures from 1000 to 20,000,K 
‘ | “ and the interval of pressures from 0.001 to 1000 atm, in (19) 
1000 4000 7000 10000 12000 
Fig.9 Molecular weight of air as a function of the temperature References 
at p = 0.001, 1 and 1000 atm. The value 4 = 28.966 is shown 1 Moore, C. E., ‘Atomic Energy Levels,” 1949. 
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(Continued from page 64) 
In a typical rocket motor (3) 


reactions. 
= 0.061 b = 3.75 
= 0.052 AH/RT? = 6.06 X 10~' for 1900K 


Using this data, dz,/dT is found to be —2.66 X 10~* mole per 
deg. The approximate solution of Equation [11] yields 
zo. Accordingly, we conclude that the reaction is frozen 
at chamber exit in agreement with the conclusion of 
Penner (8). 


An Approximate and Quicker Procedure for Estimating the 
Nature of the Flow 


This is based on an approximate solution of Equation [2] 
which enables the computation of the composition at any 
instant. We shall illustrate this procedure by considering 
the reaction 


ky 
Nz + 2NO, = Nz + N20, iF 


tecently Wegener has experimentally studied this reaction in 

a supersonic nozzle (4), and hence it is suitable for comparison 
between theory and experiment. Let c and d be the initial 
concentrations of nitrogen peroxide and nitrogen, respec- 
tively. Let y be the amount of NO, decomposed = any 
instant. Then j 


dy 


y 
At equilibrium 


Ye 


where the subscript e stands for equilibrium value. Hence 


we can write 


dy ca dt [12] 


(ye — y) (C2? — yye) 2c — 


where k, the rate of decomposition of N2Ox, is given by (5) 
ky = 1016 - ¢—13000/RT jitre mole~!-sec™! 


Squation [12] can be integrated provided y, and the cooling 
rate are assumed to be constant. There would be no serious 
error in assuming the constancy of cooling rate across the 
nozzle. y., in general, will not be constant, but an average 
value can be taken for the purpose of integration. A better 
vay would be to perform integration in successive steps, 
vhere y, may be reasonably constant. The integral on the 
right-hand side of Equation [12] is difficult to evaluate by 
rdinary methods. Graphical integration may be adopted, 
ut in the present case an approximate solution was found 
by plotting log &, against log T using the relation (5) between 
T and k, between the region of 400 and 200 K, which is in- 
volved in the experiment. 

Fig. 1 shows that for the decomposition of N,O, between 
these temperatures, the curve is a straight line to the first 
approximation. It follows, therefore 


yhere Q and mare constants. From the graph it is found that 


= 23.5 Q = 1.168 X 10-5 
Substituting the value of k, from Equation [13] into Equa- 
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= 0.0655 90) 


66 | 
54 
4 


23 24 25 26 

Log 
Fig. 1 Approximate functional relationship between rate con- 
stant and temperature 


tion rhe and integrating, we obtain 7 
Pa e e d e 
Ye—-y 2 c—Yy m+1 


[14] 


The constant of integration was evaluated from the condition 


that y = 0 when J’ = 400 K, which corresponds to chamber 
temperature of the experiment. From Wegener’s data we 
have 

c = 6.09 X 10~‘ mole 
1... d = 6.029 X mole cc 


ye = 5.48 X 104 mole/litre 


In his experiments, the cooling rate was of the order of 10® C 
per sec. Substitution of these values in Equation [14] yields 
y = y.. Since the average value of y, would lie between its 
corresponding values at 400 and 220 K, the concentration of 
NO, actually present will really be intermediate between the 
concentrations predicted by assuming near-equilibrium and 
frozen flows. This prediction agrees satisfactorily with the 
experimental results of Wegener. 

The method given above is approximate to the extent that 
an average value of z, is used. Nevertheless, it can give an 
approximate idea of the composition at any cross section of 
the nozzle. 
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Effects of Diffusion and Chemical 
Reaction on Convective Heat Transfer’ 


DANIEL E. ROSNER? 


AeroChem Research Laboratories, Princeton, N. J. 


The notion of an “effective”? thermal conductivity for 
reacting gaseous mixtures is applied to the prediction of 
convective heat transfer rates when dissociation, diffusion 
and atom recombination take place within a laminar 
boundary layer in quasi-equilibrium. The result for the 
forward stagnation region of a cool blunt nosed body is 
compared to recent solutions of the boundary layer equa- 
tions with particular regard to the explicit dependence on 
the Lewis number. 


HE THEORETICAL and experimental work in (1 


through 4)* has indicated that existing convective heat _ 


transfer data can be applied successfully to situations in which _ 


diffusion and gas phase chemical reaction take place at quasi- 3 


equilibrium provided that the ordinary thermal conductivity — 
and specific heat are replaced by suitably augmented values 
which account for these physical and chemical effects. The 
explicit consequences of an ultra-simplified version of this 
approach are briefly examined here for the case of stagnation 
point heat transfer with a partially dissociated free stream. 
The computational work of (5) is then available for com- 
parison. 

In a partially dissociated diatomic gas, if thermal diffusion 
and other secondary diffusion processes are neglected, the 
energy flux vector can be written in the ‘‘undissociated” form 


Jo = —),, grad T [1] 


where the “equilibrium” thermal conductivity A,, is related 
(2) to the ordinary “frozen” conductivity \, through 


1+ iter | (G2) - 1] 


To — this result to convective heat transfer problems, we 
note that ordinary low speed heat transfer data may be ex- 
pressed in the Stanton form . 


(Sty) Cp, ave( = Tw) [3] 


where (for gases) the dimensionless Stanton number is usually 
a power function of the Prandtl number. When dissociation, 


diffusion and atom recombination occur, the average specific 
heat and Prandtl number alone will then embody the principal 
physical and chemical effects. This naturally suggests the 
application of Equations [2 and 3], but, for the purpose of 
obtaining an explicit correlation equation, we further intro- 
duce 


(Coy ave (he hw)/(T Tw) 
It is through the ratio of these two quantities that a “chemi- 
cal” energy parameter of the form a.Q/(h.e — hw) explicitly 
enters this class of problems. Furthermore, if the degree of 
dissociation corresponding to the wall temperature is negligi- 
ble, the ratio of equilibrium to frozen Prandtl number then 
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number dependence. 


aQ 
{1 + [(Le,) 
This approach is now applied to an important particular case. 
For the laminar flow at the blunt nose of a body of revolution, 
St, can be approximated by Sibulkin’s (6) incompressible 
fluid formula, applied to the stagnation region behind the 
normal shock 


St, = [6] 


where, for heat transfer purposes, p.u. has been introduced for 
(pu)ave. Combining Equations [4, 5 and 3], we then im- 
mediately predict that the equilibrium heat transfer rate 
should be given by a correlation equation of the form 


It is now eel to compare this conclusion with recen 
numerical solutions to the complete laminar boundary layer 
equations. For the comparable case, Fay and Riddell (5 
have correlated their results with the equation 


— he 


= 0.763 — hw) X 


he 
Apart from the familiar correction factor for variable py, 
formulas [7 and 8] differ mainly in the form of the Lewis 
In Table 1, the two Lewis number 
“augmentation” factors are compared for three values of 
the “chemical” energy parameter a.Q/(he — hw) and six 
values of the Lewis number. 

With regard to the effects of Lewis number, it is seen that 
Equation [7] yields results which are higher than Equation 
[8] by less than 2.75 per cent for Lewis numbers in the neigh- 
borhood of 1.4; however, the departures increase with in- 
creasing Lewis number. Equation [7] also leads to the im- 
portant conclusion that if recombination does take place, the 


Table 1 Lewis number augmentation factors for con- | 

vective heat transfer 

a.Q 
—— Le; Eq. [7 Eq. [8 
(7] q. [8] 
1.0 1.000 1.000 
1.2 1.030 1.025 
0.25 1.4 1.059 1.048 
1.6 1.088 1.069 
1.8 1.116 1.089 
2.0 1.143 1.108 

1.059 1.050 
0.50 1.4 1.116 1.096 
1.6 1.170 1.138 
1.8 1.224 1.179 
2.0 1.275 1.217 
1.0 1.000 1.000 
12 1.088 1.075 
0.75 1.4 1.170 1.143 
— 1.6 1.250 1.208 
1.326 1.268 
2.0 1.399 1.325 
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resulting heat transfer rate, for the Lewis numbers shown, is 
about the same whether the recombination takes place at 
equilibrium within the gas phase or at a “fully” (7) catalytic 
surface. Using the result of (7) and Equation [7], this is 
seen by computing values of the ratio 


{1 + [(Les)® — 1][(@Q)/(he — hw)}} 
{1 + [(Les) — 1]l(@Q)/(he — 


This ratio never departs from unity by more than 1.4 per cent 
over the range of parameters shown in Table 1. In contrast to 
the comparable ratio 


suggested by the work of Fay and Riddell, the ratio [9] has 
the interesting property of tending toward unity for all 
values of the Lewis number when a.Q/(h.e — hw) > 1, and it 
‘s always slightly Jess than unity. 


[9] 


Nomenclature 


Cy = specific heat at constant pressure 

Di. = binary diffusion coefficient for atoms in the gas mixture 
j = sensible or ‘frozen’ enthalpy of the mixture 

= enthalpy of the mixture 

flux vector 

“frozen” Lewis number assumed constant 
Prandtl number 

heat of recombination, assumed constant 

Stanton number 


~ 
ll 


= absolute temperature ay 
= fluid velocity relative to body 

= mass fraction of atoms in mixture 

= inviscid velocity gradient at the nose 
= density of fluid - 
= thermal conductivity of mixture 

= absolute viscosity of fluid 
grad = gradient operator 


as 
| 


Subscripts 


avg = average 
= at outer edge of boundary layer 
= equilibrium 

= chemically ‘‘frozen”’ 

pertaining to energy transport 
= at the wall (surface of cool body) 
= pertaining to ordinary thermal conduction 


® 
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Uncertainty in a Calculated Specific 
Impulse Due to an Uncertainty in the 
Heat of Formation of the Propellant 


> L. J. GORDON! 
Corp., Sacramento, Calif. 


HE CALCULATED specific impulse of a propellant is 

necessarily dependent on the heat of formation of its com- 
ponents. Although these heats of formation are fairly accu- 
rately known for most propellant materials, there is a reason- 
able uncertainty in many values.? Fortunately, however, 
the uncertainties in the heats of formation of propellant 
components are small compared to heats of combustion of 
the propellants so that the errors in the calculated specific 
impulses are small. A relation between the uncertainty in 
the heat of formation and the uncertainty in the specific im- 
pulse in terms of parameters usually calculated for a particu- 
lar propellant is 


where 
I, = specific impulse a 
AH; = enthalpy of formation of propellant . 
T. = chamber temperature 
T. = exhaust temperature 
K = appropriate energy conversion factor ——. 
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and 


A derivation of the equation, a set of units for K, and a 
numerical example is given as follows. 

Although specific impulses are often calculated from the 
common approximate formula involving the heat capacity 
ratio, average molecular weight, etc., more accurate values 
are obtained from enthalpy-entropy calculations. In these 
calculations it is assumed that combustion is isenthalpic and 
the expansion process is isentropic. The path of such a proc- 
ess can be depicted on an enthalpy-entropy diagram such as 
Fig. 1, where the curves ab and cd represent the chamber 
and exhaust isobars. Assume that the uncertainty in AH, 
is such that the paths 1-2 and 3-4 in Fig. 1 represent the ex- 
pansion processes for the two enthalpies of formation AH, 
and AH, + 6AH, respectively. 

The respective heats released to become kinetic energy are 
then 


Q =H, — 


From Fig. 1 it is seen that 
6Q = (H3 — He) — (Hs — Hz) 


At any point on the H-S diagram realize that the slope of a 
constant pressure line 


(3) 
dS 
so that 
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The relationship between Q and specific impulse is 
I, = 29.5Q"/? Ibf-sec/Ibm 


when Q and AH, are given in units of keal/100 gm. There 
fore 
dl, 29.5 T. _ (29.5)? 1 


an example of the uncertainty, one finds that when 


I, = 270 lbf—sec/lbm 
T, = 3500 K 
T, = 2200 K 


an uncertainty of 20 keal/100 gm in the heat of formation of : 
component which is 40 weight per cent of the propellant 


yields an uncertainty in /, of 
“Ibm 


(29.5)? 2200 
2(270) (1 - 
High temperature shifting equilibrium heat capacities ar: 
about 80 cal/100 gm-deg K, so that the corresponding un 
certainty in 7’, is toe: 
a number small enough so as to not seriously affect the accu- 
racy of using a constant value of [1 — (7'./T.)] over the short 
6AH range. 

That the fractional uncertainty in J, is independent of the 
expansion ratio, P,/P., may be shown by rewriting the fore- 


going expression 
29.5)? 
dl, = ( 9.5) (: dAH; 


8000/80 = 


and substituting into it the approximate expression for J, 
I, = 29.5 {C,T.[1 — 


and as long as the average heat capacity C, does not var) 
with the expansion ratio, and 7’. is not changed appreciably, 
adI,/I, remains constant. 

It is also apparent that the mathematic derivation used here 
is applicable to calculating the effects on specific impulse of 
heat losses and of the propellant initial temperature. 


Transient and Steady-State Behavior 
of an Idealized Solar Powered Heat 


Exchanger 


Chance Vought Aircraft, Inc., Dallas, Texas 


The differential equation describing the transient and 
steady-state behavior of an idealized solar heat exchanger 
is developed. The equation is solved for the three possible 
operating conditions, and the steady-state temperature is 
defined. 


OLAR energy is an important source of power for space 
vehicles. To collect and utilize this energy, a typical 
power generation system would include a collector, heat ex- 
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changer and some device for converting thermalenergy into use- 
ful work. Abrupt changes in vehicle energy level require- 
ments and discontinuities in collection of the solar energy have 
made the transient and steady-state behavior of the heat ex- 
changer an important factor in the system operation. 

If the solar heat exchanger can be assumed to have a uni- 
form temperature at any instant in time and if the vehicle 
power requirements P are maintained at a constant level, the 
following expression describes the thermal behavior of the heat 
exchanger? 


Qnet = Qin (solar) out (radiation power output 


which written in differential form is 


7 dh 
Wi = SAna(1 = am) —P {1} 


2 The heat exchanger is considered alone as the reference body 
and does not include the weight of any working fluid which may 
pass through or around it. 


= 
Fig. 1 Enthalpy-entropy diag 
Solving for 6Q results in 
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We note that 


CdT 


dh = 


Since the specific heat for solids is generally a function of tem- 
perature, it is assumed that 


Ci + C3T? + C,T? 


Integrating Equation [2] for Case (c) 


21%" 


1 
- 


Wr 


] + C, [In T, — In rit [5] 


where Ci, C2, C3 and C, are determined experimentally. Hence Steady-state operation would be possible only for the con- 
Equation [1] becomes — ditions of Case (a), since both (b) and (c) result in the system 
temperature ultimately going to or approaching absolute 
= P LG i zero. In Case (a) an equilibrium state can exist, and the sys- 
ar SAnax(1— am) _ tem approaches this condition asymptotically as dT’/dr ap- 
SS és proaches zero. By setting d7’/dr equal to zero in Equation [1], 
dr Wi (Z) *) a > 3 | = the equilibrium temperature of the heat exchanger can be 
C4 + + GT ound. Hence 
hree possible modes of operation for the heat exchanger 
are indicated from Equation [2]. They are when = of SAmax(1 — Om) 1 7 
(a) Amal Qm) = S 4 a(1 Om) 
(b) SAman(1 — Gm) <P m 


(c) — an) = P 


Case (a) corresponds to transient operation which ulti- 
mately stabilizes at some equilibrium temperature 7,. Case 
b) corresponds to a transient operation, which if continued 
and were physically possible, would result in the heat ex- 
changer temperature going to absolute zero. Case (c) is 
equivalent to simple radiation cooling, where the tempera- 
ture approaches zero asymptotically. Integrating Equation 


Although the expressions above are true only for an ideal- 
ized situation, they provide insight into the general behavior 
of an actual solar heat exchanger, an understanding of the 
important parameters involved, and an approximate method 
for calculating performance. A study of this idealized behavior 
will facilitate more efficient design and provide a greater 
appreciation of the problems to be faced. 

It should be noted that the assumption of instantaneous 
temperature uniformity becomes more valid if the heat ex- 


[2] for Case (a) 
+ NT, ‘1 + 
1 1 — NT, = N?°T 
ere 
P 


~ SAmOx(1 — Om) 


Integrating Equation [2] for Case (b) 


+ MT; +1 
+ MT: V/2+ 1 2V72C: 
Ar = Ci — ] — 2— ft. 
MT; V2 MT: V5 2C., |1 + 
M? 1 — M°7;? 1 — | 1+ 
-wa 
> 
where 
a changer is small, the transient operation time is large and/or 
wh fi ‘ me - if the thermal diffusivity of the heat exchanger material is 
Nomenclature | 
1 — Gm) + ap, = solar absorptivity of solar heat exchanger 
+ am = solar absorptivity of solar collector 
An = external area of solar powered heat exchanger 
i Ax = normal projected area of solar collector 
. Cn = constants used in equation for specific heat of heat ex- 
Ti changer 
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C = specific heat of heat exchanger = power generation term : 

dT/dr = rate of change of temperature with respect to time iis = Stephan-Boltzmann constant, 0.174 X 10-® Btu/hbr 

€h = surface emissivity of solar powered heat exchanger ft? Ré 

M = general modulus used to describe behavior of solar 8 = = — constant at 9.3 X 10’ miles, 425 
heat exchanger under nonstabilizing transient , u/hr 
conditions 4 = time, hr 

N = general modulus used to describe behavior of solar Tr = temperature, R : 
heat exchanger under stabilizing transient conditions Wh = weight of heat exchanger, lb 


Steady-State Behavior of Extended 
Surfaces in Space 


J. W. TATOM! 


Chance Vought Aircraft, Inc., Dallas, Texas 


The differential equation describing the thermal be- 
havior of a uniform extended surface receiving and losing 
heat by radiation alone is developed. The equation is 
solved for three sets of boundary conditions, and a solution 
is indicated for several others. 


HE TWO heat transfer phenomena, conduction and radia- 
tion, working together establish the equilibrium tempera- 
ture of many systems in space. Hence, a study of the two in 
combination is important and necessary to the solution of 
numerous present day and future heat transfer problems. 
Consider the steady-state behavior of a rod? protruding 
from some body at temperature 7) into a solar radiation field. 
Assuming that the conduction heat flow is one-dimensional, 
the physical properties of the rod are constant and the effect 
of neighboring bodies is neglected, the following differential 
equation may be used to describe the thermal behavior of the 


rod 
{1] 


kA 
This equation can be integrated to give a 


dT 
X= pe rs 


2SaD(sin 
5kA 


kA 
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2 The following material can be applied to any constant cross- 
sectional area extended surface. The rod was used for illustrative 
purposes alone. 


aT 
dX? 


oeCT! 


kA 


[2] 


This general integral cannot be formally evaluated, «and 
numerical methods must be employed in the solution. How- 
ever, for the special case of an infinitely long rod with @ = 0 
(Case a, Table 1), Equation [2] can be integrated. The re- 
sulting solution is 


= + W9oeC/10kA 3] 


Gr-0 = V(2/5)kKAoeC (To)? 4] 


and 


This solution is applicable to rods of relatively large effoc- 
tive length (defining effective length as Le = VW BeC/KA 1) 
because their temperature distribution approaches that of an 
infinite rod. 

In Table 1, Case (b) can be evaluated easily by the substitu- 


tion ¢ = and Equation [2] becomes 


5kA 


Applying Equation [3] as a first approximation, and using the 
tabulated values of the integral presented in Table 2, 
iteration procedure can be used to determine 7. Knowing 
T,, and noting that 

dg dt 


the intermediate values of temperature can be determined. 
_ By making the substitution ae 


2kC 


Case (c) can be reduced to a form which can be integrated. 
The temperature for large values of X along an infinite rod 

receiving energy from solar radiation does not tend to zero, but 

approaches some equilibrium value. The equilibrium tem- 


: Table 1 Values of the integration constant K for @ = 0 and x/2 
Rod Boundary condition K 
(b) finite, insulated end, @ = 0 kA dXe-1 0 5 KA 
(c) finite, uninsulated end, = 0 —kA a7 = geAT’;' _ 5oeAT 
5 kA 2kC 
aX ge 5 kA Loe 
fi rsulat = _ 2 L L 
(e) finite, insulated end, 6 = r, kA 0 5 EA +. tA 
(f) finite, uninsulated end, 9 = 1/2 > | BaDT, 
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a Table 2. Tabulation of the integral 


1 

1.00 00000 1.75 0.62083 
1.01 0.08915 1.80 0.63313 
1.02 0.12616 1.85 0.64452 
1 03 0.15407 1.90 0.65512 
1 04 0.17653 1.95 0.66498 
1.05 0.19705 2.00 0.67432 
1.06 0.21399 2.10 0.69126 
1.07 0.23076 2.20 0.70614 
1.08 0.24570 2.30 0.71937 
1.09 0.26049 2.40 0.73133 
1.10 0.27351 2.50 0.74205 
1.12 0.29756 2.60 0.75168 
1.14 0.31957 2.70 0.76053 
1.16 0.33960 2.80 0.76850 
1.18 0.35773 2.90 0.77582 
1.20 0.37477 3.00 0.78255 
1.25 0.41270 3.20 0.79441 
1.30 0.44542 3.40 0.80454 
1.35 0.47399 3.60 0.81325 
1.40 0.49916 3.80 0.82084 
1.45 0.52161 4.00 0.82754 
1.50 0.54228 4.20 0.83343 
1.55 0.56058 4 40 0.83866 
1.60 0.57770 4.60 0.84332 
1.65 0.59316 4.80 0.84750 
1.70 0.60759 5.00 0.85126 


perature is governed by the expression 
SaDsin 0 = oeCTt.. [8] 
and therefore 
@ = [(SaD/oeC) -sin = U4 [9] 


Hence for a rod in the vicinity of the Earth receiving solar ra- 
diation perpendicularly and assuming a = € 


‘Tr. = 528.04R 
> 


SaD 
U= [ sin 


when substituted into Equation [2] for Cases (d), (e) and (f) 
respectively, results in the following expressions 


The 


atio 


x- 5kA aT 
= [10] 
2SaD V(1/U)Ts — 5T + 
x = 5kA dT 
Jr — 7,5) — KT — T,) 
x 
2SaD 
dT 
[12] 


V(1/U TS — — — Ty) + 50eDT,8/8kU 


Rop (@: %) 
< 
| = 
a AFinité Rdo InsucaTEO 
= 
= 1000 | 
Qa Lj: 3 Fler 
Infinite Rao 
oo} | | | 
500 | | | = 


20 
Distance Atonc ~ Fr 


Fig. 1 Temperature distribution along a 3-ft rod 


which can be solved by numerical means together with, in 
Cases (e) and (f), an iterative technique. 

To demonstrate numerically these expressions, Fig. 1 
presents the temperature distribution along a rod treated as: 
Infinite (9 = 0); finite, with an insulated end (6 = 0), and 
infinite (@ = 2/2). 

Thus the thermal behavior of simple extended surfaces in 
space can be predicted by using the previously developed 
equations. Application of the procedures and techniques to 
more complicated situations should provide a basis for exact 
calculations. 


Nomenclature | 


a = solar absorbtivity of a body 
A = area through which heat can be transferred by con- 
duction 

B = dimensional constant in definition of effective length 

C = circumference of extended surface 

dT/dx = rate of change of temperature with respect to distance 

@?T/dx? = second derivation of temperature with respect to 
distance 

D = diameter 

€ = surface emissivity 

k = thermal conductivity 

K = general constant in extended surface radiation-con- 
duction expression 

L = length 

q = rate per unit time per unit area at which heat is 
transferred 

r = constant used in integral of finite rod with uninsu- 
lated end 

= Stephan-Boltzmann constant, 0.174 Btu/hr 
ft? 

S = solar radiation constant at 9.3 < 107 miles from the 
sun, 425 Btu/hr ft? : 

U = general constant for surface in radiation field 

= temperature, R 

6 = angle between an extended surface and the sun’s rays 

X = distance 

¢ = dimensionless parameter, used in a of ex- 


| 
| 
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How 
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_ Analytical Estimates for Optimum 
Transfer Paths 

H. MUNICK' and R. McGILL? 
Grumman Aircraft Engineering Corp., Bethpage, N. Y. 


paper essentially extends the work of Lawden (1).* 
The problem considered is the minimization of characteris- 
tic velocity, or the amount of fuel consumed, using two 
impulses. 

Lawden studied the case of transfer from an elliptic orbit 
to a higher energy circular orbit. Consider a space vehicle 
traveling in an elliptic orbit with velocity components (wu, v) 
resolved perpendicular to and along the radius vector. Sup- 
pose at some point on this elliptic orbit the space vehicle is at a 
distance a from the focus, with velocity components (uw, v%). 

Suddenly an impulse is applied, giving rise to new com- 
ponents of velocity (uw, v;). The vehicle then goes into a new 
elliptical orbit, the transfer ellipse. Upon arrival at the 
higher energy circular orbit, a second impulse is applied, 
correcting the space vehicle’s velocity components to those of a 
circular orbit at a distance 6 from the focus. The characteris- 
tic velocity is given by 


V (ur — wo)? + (1 — + [we — (u/B)'2 + [1] 


Lawden then defines dimensionless parameters x, y by 


=2Vu/a n= y {2] 
This transformation reduces the problem to minimizing 
V(r — 20)? + (y — yo)? + V2? — + y? + 3r — 2 [4] (X, Yo) — (X,Y) 
subject to SS ie 
(l—r)2+y? >Al—r) [5] 


where r = a/8. 
For the case where r is approximately unity, an analytic in a RN x 


estimate can be given. An example of this is transfer from 4 
anywhere on a nearly circular orbit to a higher energy neigh- - ms) 
boring circular orbit (see Fig. 1). A second example is trans- : a Po 
fer from the near apogee region of an elliptic orbit to a higher 
energy circular orbit (see Fig. 2). Fig. 3 Optimum solution for r slightly less than unity 


Case 
Consider a small quantity e > 0. Let 


r=l-e [6] 
Rewrite Equation [4] as ee 
V(x — 20)? + (y — yo)? + 
V(2 — + + (87 — [7] 


Using Equation [6], and omitting terms in ¢ of order higher 
than one, Equations [7 and 5] become 


V(x — 20)? + (y — wo)? + Viz — 1+ (3/2) + [8] 
+ (y*/2e) 1 [9] 
Equations [8 and 9] can be given a geometric interpretation. 


The first member of [8] represents the distance from (z, y) to _ -_ a= 


Fig. 4 Transfer from outside of circular orbit to circular orbit 


(xo,yo), the second the distance from (z,y) to [1—(3/2)e, 0]. 
The constraint given by Equation [9] says (z,y) must 
lie on or outside the ellipse given by the equality in Equa- line connecting (20,yo) to [1 — (3/2) «, 0] gi ss 
tion [9]. As is presented in Fig. 3, any point on the dashed mum characteristic velocity. Once (opt, Yopt) is found (2,09 

Vi,opt) is given by Equation [2]. From Lawden’s results the 


pers ng Mathematician. Member ARS. second impulse is determined. “ 
* Systems Section, Research Dynamicist. A second case of interest occurs for r slightly greater than 


3 Numbers in parentheses indicate References at end of paper. unity (see Fig. 4). 
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Fig. 5 Transfer for r slightly greater than unity 


Case 2 J 


Consider a small quantity > 0 
r=l1+e [10] 


To a first-order approximation Equations [8 and 9] become 


Via = + (y — + 


Vix — (1+ + 
[12] 


The solution to this minimum problem, owing to constraint 
[12], is given by some (2, y) in the region (see Fig. 5) between 
the two branches of a hyperbola. There are two subcases 
here of interest: 

1 Coordinate (x,y) lies inside the branch containing 
{1 + (3/2) e, 0]. The minimum is given by the point of 
nearest tangency of an ellipse with foci at (20, yo), (1 + (3/2) €, 
0), and the branch of the hyperbola (see Fig. 6). Courant 
and Robbins (2) have treated this type of minimum problem. 

2 Coordinate (2,,yo) lies anywhere in the plane outside of 
the region given by subcase 1. For this situation the answer 
is found by joining (x,y) to [1 + (3/2) €, 0] and choosing any 
(x, y) on this line lying in the region (see Fig. 7). 

A different case of interest is for r close to zero. An example 
of this is transfer of a space vehicle near the Earth’s surface to 
a higher energy circular orbit at a considerable distance from 
the Earth’s surface (see Fig. 8). 


Case 3 


Consider a small quantity e > 0. Let 


r=€ {13} 


Neglecting terms in ¢ higher than first order, Equations [4 and 
5] become 


V — + (y — + V2? + 24 3e [14] 
z+y?—2 —2 [15] 


Lawden concludes in his analysis that the answer to the 
minimum problem occurs on the boundary given by the equa!- 
ity in Equation [15]. This reduces Equation [14] to 


V(x — x)? + (y — w)? + Ve [16] 


From Equation [16] we observe the answer is given as the 
intersection point of the normal from (29,4) to the circle given 
by the equality in Equation [15]. 

Lawden’s results are confirmed if a geometric interpreta- 
tion is given to Equations [14 and 15] (see Fig. 9). The first 
member of Equation [14] gives the distance from (x, y) to_ 
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Fig. 6 Optimum transfer for (x0, yo) lying inside branch of 


hyperbola 
(x.Y) : 


Fig. 7 Optimum solution for (x, yo) lying outside branch of 
hyperbola 


m. 


Fig. 8 Transfer from near Earth’s surface to a distant circular 
orbit 


ea 


Fiz. 9 Optimum solution for transfer from near Earth’s surface 
to a distant circular orbit _ = 
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(x,yo), the second gives the length of the tangent from (z, y) 
to the boundary curve given by the equality in Equation [15]. 
The normal distance is less than the sum of any other two 
distances. 


Acknowledgment 
namics Group for his efforts in helping complete this problem. 


Nomenclature 


u = velocity component perpendicular to the radius vector 


q 


v = velocity component along the radius vector 
a = distance from focus of point on elliptic orbit 
8 = radius of higher energy circular orbit 

uw = gravitational constant 

y 


= dimensionless parameter related to u 
= dimensionless parameter related to v 
= ratio of a to B 

small positive quantity 


e= 
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LOUIS B. WADEL' 
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The projected orbit of a 24-hr Earth satellite is a “‘figure 
8°’ oriented north-south and centered at the Equator. 
This path is calculated in terms of latitude and longitude 
as a function of the angle made by the satellite’s orbital 
plane with the Earth’s equatorial plane. If this angle is 
zero, the path reduces to a single point on the Equator. 


ppd nese in a circular orbit about the Earth (assumed 
spherically symmetrical) will have a period of rotation 
equal to the Earth’s rotational period if its altitude above the 
surface of the Earth equals approximately 22,000 miles (1).? 
If the orbit lies in the plane of the Equator and the satellite 
moves eastward, the projection of the orbit upon the Earth’s 
surface consists of a single stationary point on the Equator. 
If the orbit makes an angle a with the equatorial plane, the 
orbit’s projection is a closed “figure 8’’ path, oriented north- 
south, centered at the Equator which it crosses at 12-hr 
intervals, and reaching extreme latitudes of +a. Detailed 
properties of this path are calculated below as a function of a. 
The projected path will be obtained by first calculating the 
projected orbit’s great circle path in terms of latitude and 
longitude for a fictitious nonrotating Earth, and then modify- 
ing the results to account for the Earth’s rotation. The 
general equation for a great circle is (2) 


cos @ = sin cos + sin sin 


where 
6 = east longitude 


¢ 
A = 


Crossing of the Equator may be assumed to occur at 6 = 0, 
and the resultant path shifted later to any other crossin 

point desired. The Equation [1] must be satisfied by the 
point 6 = 0, ¢ = 90 deg, requiring c; = 0, and consequently 


[2] 


90 deg — A 
north latitude 


c tan sin = 1 
Differentiating Equation [2] and setting d¢/d@ = 0 yields 
co = cot dm = tanA, = tana [3] 


where ¢,, is minimum @ obtained, corresponding to dm, 
maximum latitude reached. The great circle Equation [2] 
can now be expressed as 


tan a tan ¢ sin 8 = 1 
Received July 1, 1959. 
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Projected Orbits of 24-hr 
Earth Satellites = tan a cot A sin 6 = 1 


On the surface of a sphere the equation for differential 
element of arc length ds is given (2) by 
(ds)?/R.? = [1 + sin? 


where R, = radius of the sphere (Earth). 

Equation [5] may be divided by (dt)? and, for the present 
case of constant speed travel, ds/dt set equal to a constant V. 
Then, utilizing the relation [4] between ¢ and 0, Equation [5] 
may be integrated to yield 


sin = sina sin (V/R,)t = sina sinw,t = sina sin 2mt [6] 
where 


w,. = angular velocity of Earth in radians per day = 27 
t = time in days measured from time of Equator crossing 


Equation [4] may be rearranged as ~ 
sin 6 = cot a tan [7] 


Once a is specified, latitude can be calculated as a function of 
time from Equation [6], and then longitude (for a nonrotating 
Earth) calculated from Equation [7]. 


The actual geographi- 
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Fig. 1 Projected orbits 
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cal path projected on the surface of the rotating Earth has 
the same latitude as calculated above, but the longitude 0, 
must be obtained from 


6. = we = — (8) 
Note also that the angle the projected path makes with the 
Equator is not a@ for the rotating Earth case. 

By differentiating Equation [8], setting d@,/dt = 0, and 
utilizing Equations [6] and [7], 6,n, the first extreme value of 
6,, is found to occur at t,, where 


sin @ 


and the corresponding value of 6,,, may be obtained from 


Equations [6,7 and 8]. Knowing that \,, = a@ (occurring for 
6, = 0) and calculating @,,, (and the corresponding ) as a by- 
product) as described previously allows quick sketches to be 
made of any projected paths of interest. Conditions of 
symmetry reduce the problem to calculating the first quarter 
period. Fig. 1 shows several such paths. The special case of 
a = 90 deg (i.e., a polar orbit) is not included by Equation 
[4] which thereupon degenerates, but may be found (for the 
first quarter period) by setting 6 = 0, so that 6, = —2rt, and 
dh = 


2rt. 
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Satellite Librations of Large Amplitude 


W. B. KLEMPERER! 


Douglas Aircraft Co., Inc., Santa Monica, Calif. 


PROBLEM of small amplitude libration oscil- 
lations of satellites or pendulous devices aboard satellites 
has been treated in Astronautica Acta (1),? and it is further 
pursued by two Technical Notes in the present issue (2). 
How such oscillations behave when their amplitudes grow 
large has been discussed by Schindler (3). The purpose of 
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the present remark is to add to the discussion the relatively 
simple exact solution of the simplified differential equation of 
motion of a dumbbell-shaped instrument of dimensions small 


-compared to the radius of a circular orbit in which it is 


carried aboard a satellite and allowed to swing in the plane 
of the orbit. 

The equation of motion then takes the form ¢ = — 30? x 
sin ¢ cos @ (where ¢ is the angular excursion from the local 
vertical and Q the angular velocity of the orbital motion). 
It is once integrated by multiplying both sides of the equation 
by 2¢ to yield the energy equation ¢? = ¢o? — 32? sin? ¢, 
where ¢y is the angular velocity during the passage through the 
vertical @¢ = 0. This immediately reveals that the max- 
imum amplitude to be attained will be given by sin @max = 
bo/V3 2 and the motion is oscillatory only if do < V3 2; 
otherwise the dumbbell pendulum will swing past the hori- 


zontal and tumble. Furthermore it is readily seen that if the 
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initial disturbance is deliberately made so that ¢@ = QQ, 
as one might want it in order to stop the absolute rotation 
for a short period while taking telescopic time exposure photo- 
graphs of distant celestial objects from aboard the satellite, 
then the swing amplitude will become ¢max = sin~! V/1/3= 
35.3 deg. This is a case for which it might also be interesting 
to compute the time required for the pendulum to return to 
¢ = O where the speed would be brought back to ¢ = 0; this 
time turns out to be equal to the time it requires the satel- 
lite to travel over two radians of its orbit. 

The solution of the time history of the pendulous motion 
is the integral of the reciprocal of the square root of the 
energy equation spelled out previously 


t= do/W do? — 30 sin? 
For ¢o = kV/3Q, when k> 1, this reduces itself to the com- 
plete elliptic integral uae 


and the result is a continuous tumbling at a cyclically varying 
rate. If k < 1, the substitution sin ¢ = & sin yw transforms 


the expression into the incomplete elliptic integral ip te 
t= (1/V3)S 


for which the solution is periodic. 

Three cases are of special interest: 

1 k= 1/3 =0.577, which means = Q and for which 
the swing to maximum amplitude takes the same time as 
for the satellite to travel one radian. 

2 k = 0.965 for which the period of the oscillation is the 
same as that of the orbit. 

3 The borderline case for k = 1 or do = 32 forwhichthe 


motion aperiodically and asymptotically approaches the 
horizontal. In this case the elliptic integral degenerates 
into the simple form 


t = (1/0/32) sec db = (1/03 Q) In tan (4/4 + $/2) 


Fig. 1 depicts the time history of the motion plotted in terms 
of excursion @ against the angle 6 = Q ¢ subtended by the 
satellite orbital motion in the meantime, for the three special 
cases just described and for several other example cases, one 
of lesser and one of intermediate amplitude and one tumbling. 

Incidentally, the motion of the dumbbell, as here described, 
is of a character closely related to that of a conventional 
pendulum of length 7 swinging from an Earth-fixed fulcrum, 
for which the equation of motion to large excursion is ¢ = 
—(g/l)sin $; this is solved by the same kinds of elliptic 
integrals, but with the argument ¢/2 instead of ¢. 

These derivations, interesting as they are, are limited to 
the dumbbell configuration. If analogous reasoning were 
applied to the case of a more compact body then the motion 
would appear slowed down by an appropriate factor. Further- 
more, the findings apply strictly to a satellite in circular or! vit 
around a central body endowed with a spherically symmeti- 
cal gravity field. If the orbit is elliptical or if the field has 
quadrupole character or possesses anomalies, pertubations 
enter into the picture. 


1 Klemperer, W. B. und Baker, R. M. L., Jr., ‘Satellite Libra- 
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Aside from an idealized circular orbit assumed by Klem- 
perer and Baker, satellite libration phenomena, and hence 
passive stability also, occur with bodies traveling on orbits 
of small but finite eccentricity. In such a case, the motion 
of the satellite about a position of stable equilibrium is 
found to be a solution to Mathieu’s equation with an 
amplitude in radians less than or equal to the orbital 
eccentricity. 


MPLICIT in the analysis of a paper by Klemperer and 
Baker (1),3 which dealt with the libration in the orbit 
plane of an idealized dumbbell and a solid ellipsoid, was a 
constant angular rotation of the dumbbell system at a velocity 
of n radians per unit time, (i.e., n stands for the mean daily 
motion). Of course, such a constant rotational velocity did 
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not affect the dynamical equations of the system, since the 
orbit was perfectly circular and the geocenter effectively 
rotated about the dumbbell exactly in step with the dumb- 
bell’s constant rotation. In the case of an elliptical orbit the 
geocenter does not keep pace with this constant angular veloc- 
ity but rather leads it from the perigee to apogee and lags it 
from apogee to perigee. See Fig. 1. This angle of lead or 
lag is simply the difference between the true anomaly v and 
the mean anomaly M of the satellite; i.e., the constantly 
rotating axis of the dumbbell makes one revolution per 
orbital period, and therefore its angle with respect to an iner- 
tial direction in space is simply M = My — n (t — t), where 
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M, is the value of the mean anomaly at the epoch &, and n 
is the mean daily motion in radians per unit time. 

As a result of their different distances from the geocenter, 
the two spheres suffer different attractions 


where 
Fg = gravitational force ; 
k, = geocentric gravitational constant 
= 
m = mass of a dumbbell sphere ~ 
r = radial distance from the geocenter 


Subscripts 1 and 2 refer to dumbbell spheres 1 and 2. 
The respective leverages a; and a2 of these attractive forces 
are now 


=([rLsin(@ +v — M)]/n [2a] 
a. = [r Lsin(@ +v — M)]/re [2b] 
as can readily be seen from Fig. 2, where L is the half-length 
of the dumbbell connecting rod and ¢ is the angular departure 
from the constantly rotating reference axis (¢ is measured 


negatively). 
Hence, the net torque due to the different attraction forces 


is 
T = Fag — Fea = 
k?mrLsin(¢ +v — M) (1/r3 — 1/r23) [3] 

The centrifugal force torques cancel out and contribute 
nothing. 

So long as the dimension L is very small compared to any 
orbital radius r, the individual radii can be approximated by 
mn =r—Lecos(¢+v— M) andr, =r+Lcos(¢+v — 
M). Introducing these into Equation [3] and expanding the 
reciprocal cubes into a series of power terms, we find that the 
torque equation becomes 

T = (6k.2mL2/r*) sin (@ + v — +v—M) [4] 

For an elliptical orbit of small eccentricity, e, r can be ex- 

pressed approximately as 
r —ecos M) [5] 
Expanding 1/r* for small eccentricity e, we find that 


= (6k,2/a?) mL? (1 + 3e cos M) sin(@ +v — M) X 
cos (@ +v—M) [6a] 


Furthermore, assuming the angle (@ + v — 1) to be small, 
we have 


= (6k,2/a*) mL? (1 + 3e cos M) (6 + v — M) [6b] 
From (2), p. 171, Equation [64] 


v — M =2esin M 


P? = or k,?/a® = (27/P)? = n? 
(P is the period of the orbit.) Therefore 
T = 6n? mL? [(1 + 3e cos M) (@ + 2e sin M)] [6c] 


where terms of order e? are neglected. 

If the dumbbell is allowed to swing in response to this 
torque, the result will be an angular acceleration ¢ (in the 
uniformly rotating reference coordinate system) governed by 


T=-2mL'?¢ [7] 


the minus sign arising because ¢@ has been counted nega- 
tively. Equating [6c and 7], we find that the part 2 m L? 
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+ 3n? K*(1 cos nt) = 


> 


46 


=L sin (6 +v-M) 
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Fig. 2 Dumbbell librations of an eccentric orbit 


cancels out and 
= —3n*(1 + 3ecosM)¢ +2esinM [Sa] 

or 
+ 3 n* (1 + 3e cos nt) = — 6n*esin nt [8b] 
If z, a and q are defined as 


M+t+n(t—h) 


+12 
(9/4) e 


2 
ds 


(My = nb) 


ll> & 


q 


then Equation [8] becomes 


which is the canonical form of Mathieu’s equation with a 
“Tf forcing function of — (32/3) q sin 2 2; see (3). 
Pad The prolate spheroid analysis follows along very similar 


Also, it should be noted that 7 pe, 


lines. Utilizing the relations from (1), it is readily seen that 
6 n? e sin nt [10] 
where 


K? = &/(2 — 

and € is the eccentricity of the prolate spheroid. Thus 

q = (9/4) i 

The solution of Mathieu’s equation follows along two 
somewhat different lines. Its application to physics involves 


finding the appropriate constant values of a (depending of 
course on qg) such that the equation possesses a periodic solu- 
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tion. Problems of this type arise in connection with tidal 
waves in a cylindrical vessel having an elliptic boundary, 
certain forms of steady vortex motion in an elliptic cylinder, 
the decay of magnetic force in a metal cylinder, etc. (4). For 
astronomical problems (such as determining the motion of 
the lunar perigee) a and g are known constants and one at- 
tempts to solve for the (usually) aperiodic function.® 

In our general discussion of libration, both of these problems 
are relevant. With the geometry and orbit of the moon or of 
the satellite specified, we can only determine the (usually) 
aperiodic motion. However, if the orbit parameters and the 
moment of inertia of the satellite can be varied then we may 
choose to seek out the constants of the system leading to 
a periodic solution. Because of the complications arising from 
the forcing functions, an analytical solution to the problem is 
not very promising and recourse to an electronic computer is 


5 Actually Hill’s equation is used in lunar theory. It is a 
slightly more general relation involving an infinite series in 
cos 2nz. See (5). 


indicated. Such work has been carried out by Stocker and 
Vachino (6) and by Robinson (7). 

Tne author wishes to acknowledge the many helpful sugges- 
tions and comments made by Drs. W. B. Klemperer, Eric 
Durand and L. G. Walters during the preparation of this 
note. 
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oo Libration of a Prolate 


Ellipsoidal Shell! 


ROBERT M. L. BAKER Jr.’ 


University of California, Los Angeles, Calif. 


A satellite in the form of a thin prolate ellipsoidal shell 
contained between two similar ellipsoids is assumed to 
resolve in a steady circular orbit around the Earth. It 
rotates essentially so as to show the same pole of its major 
axis to the Earth, except for a small angular oscillation 
(libration) superimposed on this rotation and, conse- 
quently, behaves in a similar fashion to the solid ellipsoid 
analyzed by Klemperer and Baker. 


ET A Cartesian coordinate system wow be defined as hav- 
ing its origin in the mass center of the satellite, its w axis 

coincident with the radius vector from the geocenter, while 
the v and u axes are in and normal to the plane in which an 
oscillation of the long figure axis (a — a’) of the ellipsoid 
against the radius vector axis w lies. (This plane may be 
conceived as coincident with the plane of the orbit, but the 
solution is not confined to this restriction.) See Fig. 1. At 
the instant ¢ the angular excursion between a and w shall be 
denoted by ¢. In order to break down the mass of the thin 
ellipsoidal shell into infinitesimal particles it is convenient 
also to introduce a cylindrical coordinate system X, y, 6, 
solidly and symmetrically anchored in the ellipsoid, where X 
is the axial coordinate, y the radius in the body of revolution 
which assumes the value Y; at the inner surface and Y2 at the 
outer surface, and @ is the azimuth angle counted from the u 
direction, as illustrated in Fig. 2. 

An increment of mass in the prolate ellipsoid having a 
Received July 1, 1959. ;. 

1 This paper presents the results of a study initiated at Douglas" 
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dm = [1] 
ara 


The distance of this increment above (or below) the w 
plane is 


w= Xcosd+ ysin 6 sind [2] 


See Fig. 3. 
For every small ¢ this is, with close approximation 
w= X + (ysin [3] 
Since the prolate ellipsoid is symmetrical about the vw 
plane, it suffices to consider the moment of 2dm integrated 
from 6 = — 7/2 to +7/2, and acting on the arm L which is 
the distance from the projection of the increment dm on the 
vw plane to the center of mass. If dF is the incremental force 
acting on 2dm and a is the angle of Z with respect to the v 


CENTER 
OF MASS 


Fig. 1 Coordinate system 
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Fig. 2 Symmetrical masses 


axis, then the incremental] moment = L dF cos a, where the 
leverage is, as readily gleaned from Fig. 4 


Leos a = X sing — ysin 6 cos [4a] 
which for small ¢ is approximately 
Leos a = Xo — ysin [4b] 


Analogous to Equation [2] of Klemperer’s and Baker’s 
paper,* the incremental torque due to the symmetrical pair of 
increments of mass above and below the w plane as per 
Equation [4a] (situated symmetrically with respect to the 
center of mass of the ellipsoid) is i i 


dT = 2dmrL cos a(1/r,3 — 


a” 
where 
4 > Vv 
r, = — (L/r) sin a] 
 Fig.4 Lever arm 
and 
re = r{1 + (L/r) sin a] 
Apyeuulunately The torque equation becomes 
a Y2 
1/ri3 (1/r?)[{1 + (3/r)L sin a) Id = —12n%p 28 dx dy (X?¢ — 
and Xysin 0 — y? sin? A)y 
1/r23 & (1/r*)[1 — (3/r)L sin a] [7] 


Integrating 6 from —7/2 to 1/2, we find that the sin @ 


Noting that for the circular orbit k,?/r? = n?, we find that the drops out and 7 
incremental restoring moment caused by the antisymmetrical . ¥ 7 
pair of symmetrical 2dm acting on the leverage [4a] now reads I¢ = —12n*prd f, dX f a ‘ ydy(X? — y?/2) 

1 


dT = 2dmn?-6(X cos @ + y sin sin 


(X sin — ysin@ cos) [8a] — — 
= 
dT = 12dmn? (X + yg sin — ysin 8) 

= 12dmn? (X%p — Xy sin 0 — y%@ sin? 0) [Sb] For the generatrix ellipse Y,2/b2 + X?/a:? = 1, where q 
and 6, are the semimajor and semiminor axes of the inner sur- 

* Klemperer, W. B. and Baker, R. M. L., Jr., “Satellite Libra; face. Therefore a 

tions,’’ Astronautica Acta, vol. 3, fasc. 1, 1957, pp. 16-27, an tm 4 
vol. 3, fase. 3. = — 1+2 [10 
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The integral of Equation [9] is for either surface (4/3)rp 3pab? [14b] 


(b?/a® + b4/4a*)- = = — 6%) [ll] 


{12a] 


Let the semi-axes of the inner and outer ellipsoidal surfaces be 


a, = a, bk = band a, = a(1 + p), be = W(1 + pw) 


therefore 
= —nrp(4-5/5)p aba? — b*)p [12b] 
where y? has been neglected. a 
For a thin ellipsoidal shell 
3) 


4m + b?)/3 (a? + b?) 


Ww hich i is less than the frequency of the dumbbell by the frac- 
tion V (a? — b*)/(a? + b?), the same result as obtained for the 
solid homogeneous prolate ellipsoid found in Klemperer’s and 
Baker’s paper® Equation [20]. 

Introducing the numerical eccentricity ¢ of the prolate- 
ellipsoidal shell by the notation 


A (a? — b*)/a? 


we have 


la? — bt 


| 


[16] 


Consequently, the frequency of the shell libration is given by 
w = nV/3K 


The suggestions and contributions to this note made by 
W. B. Klemperer are gratefully acknowledged. 


[17] 


Effect of Thrust Misalignment on Gyro- 


Stabilized Vehicles 
RUSSELL P. NAGORSKI? 
Consideration of thrust misalignment is a well-known 
factor in many flight control system dynamic analyses. 
This same undesired influence, if uncompensated, has the 
effect of causing the deviation of the velocity vector from 
the prescribed flight path with its attendant lateral excur- 
sion. The purpose of this paper is to present an analytical 


development of the significant expressions that effect the 
flight profile due to thrust misalignment. 


Aerojet-General Corp., Azusa, Calif. 


Conventional Gyro-Stabilized Flight Control System 


TYPICAL control system can be represented as shown 
in Fig. 1 where the component time constants, small rel- 
ative to the vehicle control period, are neglected. 

The angular displacement @ is indicated by a displacement 
gyro; the angular rate 6 by a rate gyro; the summation of 
both in turn defines the orientation of the thrust vector y. 
The orientation causes an angular acceleration of the vehicle 
to close the loop. From Fig. 1 a set of differential equations 


defining the motion can be established 
+ Kb = 7 


An expression for the vehicular axis orientation in the time 
domain can be deduced by equating Equations [1 and er thus 


A(t) = E + et bt — cos ] 


where 


Received July 13, 1959. 
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Similarly an expression for the thrust vector orientation can 
be given by 


The equation defining the time history of the lateral veloc- 
ity vector component is derived by the basic equation 


VAt) = r fw dt (5 | 


Neglecting the transient terms of Equation [4] and assum- 
ing a time-based function of mass of the form 


= e{1 — e*[cos bt — (a/b) sin bt)} (4) 


the integrated solution to Equation [5] is seen to be 
V(t) = (—Te/Moa) In (1 — at) (7| 


Subsequently, a further integration defines the lateral ex- 
cursion, which is simply 


X(t) = (Te/Mya?)((1 — at) In (1 — at) + at] [8] 


- 


Ke 


Conventional attitude stabilization system schematic 


Fig. 1 
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In order to illustrate the utility of Equations [7 and 8], the 
following typical illustration is presented using the perform- 
ance data from a typical booster stage using a 4.5ks-115000 
solid in a particular multistaged vehicle: 

Initial mass of combined stages, My) = 240 slugs. 

Vehicle mass at burnout of booster stage, 1/,, = 122 slugs. 

Thrust (assumed constant during 4.5 sec burning period), 

T = 115,000 lb. 
Possible thrust misalignment, € 
Burning time, 7 = 4.5 see. 


= 1/2deg = eee radians. 


Velocity at burnout, V = 2,300 fps. 
Therefore 
1 — at = (My — Mw)/Mo = 0.493 
and 
also 


= 0.113 


Substituting these values in Equation [7] results in V, = 
26.1 fps. Therefore the deviation of the velocity vector is 
seen to be Y = 26.1/2300 = 0.0135 radian or 0.77 deg, a 
significantly large figure. 

A similar substitution in Equation [8] would provide the 
lateral excursion with the parameters given. This is, how- 
ever, usually considered of secondary importance, since many 
of the present-day missions have relatively large coast periods; 
consequently, the trajectory accuracy is essentially dependent 
on the velocity vector attitude accuracy as well as its — 
tude. 


Modified Gyro-Stabilized Flight Control System 


If the dynamics of the system previously described were 
changed so as to include a linear acceleration feedback to the 
summing amplifier by means of a body mounted accelerom- 
eter, the schematic form would appear as shown in Fig. 2. 

By the method of derivation previously demonstrated, a 
series of equations denoting the vehicle characteristics in the 
time domain can be evolved. These are 


= + et (: sin di — cos w)| (9) 
1 


where 
c = — K,) 
—K,) 4K;%(1 — K,)? 
and 


€(1 — K,){1 — dt — (c/d) sin dt}} [10] 


Ke 

‘Fig. 2. Attitude stabilization system schematic with lateral 

Neglecting transient terms as before + . 

= [—Te(1 — Ks)/Mya] In (1 — at) [11] 


= [Te(1 — {(1 — at) In (1 — at) + at} [12] 


Here the usefulness of an acceleration feedback is illustrated, 
namely, that any finite positive value of accelerometer gain 
will tend to attenuate the undesired effect caused by thrust 
misalignment. In practice, the objective of the system en- 
-ginee r is to negate the effect of thrust misalignment by intro- 

-duci ing an acceleration feedback as shown here with a unity 
gain value. 


Nomenclature 

I = moment of inertia, slug-ft? 
K, = displacement gain constant 7 
K, = rate gain constant, sec . 
K; = dynamic constant = I/TI, sec? 

K, = acceleration gain constant 

l = moment arm, ft 

M_ = instantaneous mass, slugs 
M, = mass at time? = 0, slugs _ 
My = mass at burnout, slugs 

t = time, sec 

V = velocity along the principle axis, fps 


Vz = velocity normal to the principle axis, fps 
= normalized flow rate, 

y = orientation of the effective thrust vector, radian a 
6 = angular attitude displacement, radian 
6 = angular rate, radian/sec 

¥ = velocity vector deviation due to thrust misalignment, 


radian 
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Cutoff Impulse or ‘een Time 

N A RECENT note, Kelly (1)? discussed the effect of the 

thrust termination process on the range dispersion of a 
ballistic missile. He concluded that variations of the cutoff 
impulse and of the missile mass at cutoff both contributed to 
appreciable dispersions. However, as is shown below, it is 
comparatively simple to arrange the guidance system so that 
errors due to the first source of variability are reduced and 
these due to the second source completely eliminated. 


Suppose that the cutoff impulse, instead of being quoted 
directly as the actual nae 


is specified by the equivalent cutoff time ¢,. defined by the 
equation 


S. L. BRAGG! 
Rolls-Royce Ltd., Derby, England 


where Pi is the thrust just before cutoff. This equivalent 
time is, as it were, an overrun time. It is the time for which 
the motor would have to run at full thrust after the cutoff 
signal was given, in order to produce an impulse equal to the 
actual integrated cutoff value. 

In conventional guidance systems the nominal cutoff im- 
pulse and nominal missile mass at cutoff are used to calcu- 
late the nominal increment in final velocity (AU;)nominai = 
(Ieo)n/(Meo),. (The effect of the change in missile mass 
during cutoff is negligible.) The guidance system is ar- 
ranged to signal cutoff when the missile velocity reaches a 
value (AU,)nominai less than the required one. Deviations 
of either J., or M.. from the nominal values produce propor- 
tional error in AU,;. 

Suppose on the other hand that the guidance system is 
arranged to measure the momentary missile acceleration a 


Received July 20, 1959. 
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(as i is inevitable i in all smertiel sy stems) and is arranged to siy- 
nal cutoff when the missile velocity reaches a value at.. less 
than the required one. This value will vary slightly froin 
flight to flight, compensating automatically for the variatiors 
of missile mass at cutoff. 

For example a higher than nominal mass at cutoff will 
yield a low acceleration at that time. On the proposed sy-- 
tem cutoff will then automatically be signaled at a slightiy 
higher missile velocity than normal, so the final velocity wi! 
be correct. 

Although errors due to variations in missile mass at cutof 
are eliminated in this way, errors due to the variability of 
will remain. These, however, are likely to be smaller than 
the errors arising from the variability of J... This is because 
variability of cutoff impulse arises from two main sources 
variation of the sequencing and valve operating times, and 
variation of the thrust level of the engine at cutoff. The 
first source produces identical variabilities in Z.. and t.., but 
since the second source would be expected to increase J. in 
the same proportion as F,., it has no effect on to. 

The proposed system has also a secondary advantage. It 
is unnecessary to alter the constants involved in computing 
the correct time to give the cutoff signal when an alternative 
range (that is, missile mass at cutoff) is required. 

It is interesting to note that if we apply the analysis of 
Rodean (2) for instantaneous stoppage of propellants to « 
fixed geometry thrust chamber in vacuo, we find that 


teo = 2m./(k + 1)m 


where 
k = specific heat ratio 
Mm, = mass of gas in the chamber under steady running condi- 
tions 


3. 


corresponding outflow rate 


However, it appears unlikely that main propellant valves can 
in practice be closed in a time appreciably smaller than the 
gas residence time so this particular relation is of academic 
interest. 
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Note on “Interplanetary Trajectories 
Under Low Thrust Radial Acceleration” 


H. K. KARRENBERG! 


Systems Corporation of America, Los Angeles, Calif. 


HIS note presents corrections and additions to Cope- 

land’s? derivations of @ and 7 for the specific cases, a = 
1/9 and a = 1/6, and for the generalized case, a > 1/8. 

For the case a = 1/9, Equations [21 and 22] should re- 
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T= 20: 


3 tan 
2 of t- (1/4) sin? 6 


For a > 1/8, Equations [27 and 28] should respectivel) 
read 


spectively read 
a - 


F (6 3) + are tan 


4 1 + cos 
V1 — sin? @ — F(0, k) + 


sin 6 


+ constant 


E(8, | 
arc sin (k sin 6) — (1/V 8a)F(0, k) + constant 
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In order to conform with the previous cases the following two 
equations should be added for a = 1/6 


sin 0 


2 9 
3 3 


In the middle of the right-hand column of p. 270,? the equa- 
tion for r which indicates when escape velocity is reached 
should read 


=1+1/2a 


Nomenclature 


¢@ = angle subtended at sun between vehicle and point of de- 
parture from Earth’s orbit 
= nondimensional time 


T 
a = nondimensional vehicle acceleration 
F = elliptic integral of the first kind . 4 =| 
E = elliptic integral of the second kind af a 


Use of Inert Gases in Flame Stabilization 


P. ROY CHOUDHURY! and ALI BULENT CAMBEL? 


Northwestern University, Evanston, III. 


| IFFERENT aspects of the phenomena associated with 
the stabilization of premixed flames by gaseous jets have 
been reported previously in this Journal (1)? and elsewhere (2). 
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In earlier publications we had suggested that a small critical 
zone governs the mechanism of stabilization. and that inert 
gases may not be useful as stabilizing media. Since then we 
have performed some experiments which indicate the need 
for modifying these descriptions. 

Our studies indicate that in addition to the critical zone, 
there must be considered the influence of a diffusion and 
mixing interface as suggested in (2). Furthermore, we find 
that argon, nitrogen, carbon dioxide and helium may be used 
as stabilizer gases. However, we noticed that the behavior 
of carbon dioxide and helium was erratic. 

The data obtained with argon and nitrogen are plotted in 
Fig. 1 and compared with those for air. As may be seen, 
two different pressures and three grid sizes are compared. 
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Fig. 1 Performance of inert gases as flame stabilizers 
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Neither factor influences the performance systematically. 
The ineffectiveness of the grid size variation is understand- 
able because hot wire anemometry indicated very high levels 
of turbulence in the stagnation zone. 

The fact that inert gases may be used as stabilizers suggests 
that in addition to the chemistry, the fluid mechanics and 
the associated transport phenomena must be considered 
sriously. Probably among the principal difficulties in doing 


this are the many little known aspects of jet interaction and 
the meager knowledge concerning the turbulent and convec- 
tive terms of transport properties. 
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Addendum to “Thermal Protection of = 


Re-Entry Satellite” 


General Electric Co., Philadelphia, Pa. 


N AN EARLIER study (1),* the author presented an 

analysis of the total mass transfer required at the forward 
stagnation point of a re-entry satellite. Included were 
transpiration cooling, vaporization cooling (T = 0.25) and 
sublimation. 

In a more recent study (2), the author presented analyses of 
the ablation of quartz-like refractories, thermosetting rein- 
forced plastics, a thermoplastic material (Teflon), and graph- 
ite, under the conditions of hypersonic re-entry. 

It is interesting therefore to compare the total mass loss 
requirements for all of the cases treated in (1 and 2). 
done in Fig. 1. 

It is seen that Teflon produces the greatest, and graphite 
the least mass loss, and all of the other materials treated fall 
between these upper and lower limits. 

However, it must be noted that Teflon ablates in a quasi- 
steady fashion during re-entry and hence requires a relatively 
thin insulation layer behind the ablating layer. Graphite 
behaves primarily as a heat sink material, which also under- 
goes surface combustion, and thus the curve marked graphite 
represents only the mass ablated, not the total required for 
insulation. 

The curve marked ‘vaporization cooling [ = 0.25,” in 
(1), Fig. 5, coincides with the curve marked quartz, and 
hence is not shown here in Fig. 1. It is also remarked that 
the curve marked “transpiration cooling” in (1), Fig. 5, 
inadvertently included the weight of the porous wall while 
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1 This work was performed under the auspices of the USAF 
Ballistic Missile Division, Contract no. AF 04(647)-269. 
2? Manager, High Altitude Aerodynamics, Missile and Space 
Vehicle Department. Member A 
3’ Numbers in parentheses indicate References at end of paper. 
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ig. 1 Total mass transfer at stagnation point of re-entry 
satellite 
in Fig. 1 only the mass of air actually injected into the bound- 
ary layer isshown. This mass requirement is truly optimistic 
because of the weight of the porous wall, the storage tanks 
and flow regulating equipment required. 

Note also that although the quartz and the reinforced plastic 
seem equally competitive on a total mass loss basis, the 
depth of thermal penetration will be larger for quartz than 
for typical reinforced plastics, since the former has a highe r 


thermal diffusivity (3). 
’ ARS 
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Comment on “Adiabatic Wall Tem- 
perature Due to Mass Transfer Cooling 
With a Combustible Gas” 


D. B. SPALDING! 


N A COMMENT (1)? on a technical note by Sutton (2) 
I omitted a factor of 2 which is used by that author in 
applying the flat plate boundary layer solutions to the axi- 
symmetrical stagnation point. This error is corrected in 


ae) Imperial College, London, England 
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Fig. 1, which corresponds to Fig. 2 of (1), and which relates 
the dimensionless mass transfer rate at the stagnation point 
to the dimensionless driving force. 

Sutton’s curves, which are valid for constant properties, 
are those marked (2) and (3); the former, valid for Pr = 1, 
corresponds to Fig. 2 of (2); the latter embodies a correction 
for Pr = 0.7 (8). Comparison with curve (1), reeommended 
in (1) as the exact solution for the constant-property laminar 
axisymmetric stagnation point, shows that the error asso 
ciated with the use of the flat plate solutions is considerablh 
less than was suggested in (1). 

In plotting Fig. 1, the opportunity has been taken to intro- 
duce data presented by other authors, with particular refer- 
ence to the influence of property variations in the boundary 
layer. Accordingly the viscosity and density appearing in 
the ordinate expression require definition ; in Fig. 1 the values 
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at the outer “edge’’ of the boundary layer are employed, 
namely p, and pe. 


Curve (la) represents constant-property solutions © 
the paper by Reshotko and Cohen (4). Since these wereob- A 
tained by interpolation in the same exact solutions as those = 


underlying curve (1), the two curves naturally exhibit only 
minor differences. Reshotko and Cohen also obtained exact 
solutions of the boundary layer equations valid for variable 
properties; their results for 7.,/T.. = 0.25 are plotted as curve 
(4) on Fig. 1. 

Also plotted, as individual points, are the variable-property 
solutions of Howe and Mersman (5), who assumed: paT®’, = 
kaT®*5, c,aT°8, paT-!. The remarkable feature of these 3 
results is that, though 7.,/7., varies between 0.25 and 4, the 
points all lie close to the constant-property curves (1) and 
(la). This behavior contrasts with that exhibited by curve 
(4). The difference may perhaps be attributable to the 
different property variations assumed by Reshotko and | 
Cohen, e.g., wa7'; their assumptions are rather less in accord- 
ance with the properties of air than are those of Howe and he 
Mersman. 

The small influence of the temperature ratio on the mass 
transfer relation, when the properties are inserted at their 
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“outer edge’’ values, is also exhibited by the solutions of B 

Fay and Riddell (6). These authors considered heat trans- Fig. 1 Curves for prediction of mass transfer rate at laminar 
fer in the absence of mass transfer, a case which represents axisymmetric stagnation point. All curves are calculated for 
the B + 0 limit of the present problem. Their results may be Pry = 0.7 and T./T, = 1 unless otherwise stated 


expressed as 


= 0.957 {| —— 1 


This asymptotic relation is drawn on Fig. 1 as a broken line 

marked (fore which The 4, gn stream a ote “eg” of andy 

: ’ x = distance along surtace from stagnation point 

nonexistent if one makes the Reshotko-Cohen assumption: Pr = Prandtl number 
up is independent of temperature. T = absolute temperature of gas 

Considering the results represented on Fig. 1 as a whole, un = dynamic viscosity of gas 
it appears that the relations represented by curves (1) or (1a) p = density of gas aa 
may be used with some confidence for variable-property ; . 
boundary layers provided that u and p carry the suffix e. ieateihatea —_ 

In conclusion I would like to mention that Fig. 1 (the h-f e€ = ingasstream just outside boundary layer — 
diagram) of (1) was inaccurately plotted. A more careful w = ingas immediately in contact with surface — 
working out of the numerical example discussed in (1) has eo ae 
shown that B = 1.725 with combustion and 1.42 without. ones 

1 Spalding, D. B., ‘‘Adiabatic Wall Temperature Due to Mass Transfer 


Cooling With a Combustible Gas,’’ ARS JourNat, vol. 29, no. 9, Sept. 
1959, pp. 666-668. 


Nomenclature \< : 2 Sutton, G. W., ‘‘Adiabatic Wall Temperature Due to Mass Transfer 
4 a tit Cooling With a Combustible Gas,’”” ARS Journat, vol. 29, no. 2, Feb. 1959 
B = “transfer number” or driving force for mass transfer pp. 136-137. 
(N.B. the name “enthalpy ratio” for this quantity is 3 
ins iently general: , rittten. = 4 Reshotko, E. and Cohen, C. B., “Heat Transfer at the Forward Stag 
nesses ntly general; for ad may be written in ae nation Point of Blunt Bodies,” NACA TN 3513, July 1955. 
of concentrations and involve neither enthalpy nor 5 Howe, J. T. and Mersman, W. A., ‘Solutions of the Laminar Com- 


pressible Boundary Layer Equations With Transpiration Which Are Appli- 


temperature ) 
cable to the Stagnation Regions of Axi-symmetric Blunt Bodies,’"”, NASA TN 


Cp = specific heat of gas at constant pressure 

k ee a thermal conductivity of gas , a 6 Fay, J. A. and Riddell, F. R., “Theory of Stagnation Point Heat 
m” = mass transfer rate per unit area (e.g., in lbm/ft?h) =. Transfer in Dissociated Air,” J. Aeron. Sci., vol. 25, 1958, p. 73. 
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Jet propelled helicopter rotor (2,894,589). 
F. F. Ehrich, Philadelphia, Pa., assignor to 
Westinghouse Electric Corp. 

Ramjet or turbojet engine mounted in 
the rotor hub. Blade tips have openings 


for receiving the incoming air and for 
discharging motive exhaust gases. The 
air intake receives the supply of air for 
the engine by “ram effect’’ with attendant 
aerodynamic advantages. 


Jet propulsion device for operation through 
fluid medium (2,891,381). F. Zwicky, 
Pasadena, Calif., assignor to Aerojet- 
General Corp. 


Reaction engine with an automatically 
operable valve dividing an elongated pas- 
sageway into two parts, and located be- 
tween the inlet opening and exhaust nozzle. 
A shorter duct located in the passageway 
comprises a vestibule and combustion 
chamber. 

Guided missile (2,892,600). W.K. Ergen, 
Oak Ridge, Tenn., assignor to the U. S. 
Army. 


Receiver and Antenna Servo Drive 


Radar Transmitter 


Means cooperating with a control gyro- 
scope whereby a missile may be turned 
from its vertical path into the horizontal, 
at a predetermined altitude, without the 
application of precessing forces to the 
gyroscope, and without change of position 
of the gyro spin axis relative to the vertical. 


Epitor’s Nore: Patents listed above 
were selected from the Official Gazette of 
the U.S. Patent Office. Printed copies 
of patents may be obtained from the 
Commissioner of Patents, Washington 25, 
D. C., at a cost of 25 cents each; design 
patents, 10 cents. 
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Method for generating impulses in electric 
detonators (2,889,776). U. Gunther, Zu- 
rich, Switzerland, assignor to Inventa 
A.G. Fuer Forschung und Patentverwer- 
tung. 

Production of high induction potential 

for fuses by the interruption of a magnetic 
circuit. One leg of a U-shaped permanent 
magnet contains an explosive charge which 
is detonated upon impact with a target. 
Destruction of the magnet interrupts the 
circuit. 
Method of operating rotary wing aircraft 
including jet driven rotors (2,889,887). 
P. H. Stanley, Glenside, Pa., assignor to 
Autogiro Co. of America. 

Rotors driven at blade tip speed of 

between 680 and 760 fps controlling the 
craft by controlling the rotor blade path. 
Variations of attack angle from all causes 
is kept within the bucket range of the 
rotor blade section. 
Target assignment for radar tracking 
apparatus (2,891,244). J. J. Pastoriza, 
Boston, Mass., assignor to the U. S. Air 
Force. 

Aprlication to a tracking means of a 
pair of voltages proportional to the coordi- 
nates of the selected target and a plan 
position indicator for displaying the tar- 
gets. The assignment system consists 
of an electrically resistive overlay for the 


George F. McLaughlin, Contributor 


plan position indicator. A manually 
controlled switch applies the voltages to 
the selected tracking means. 
Signal tracking device (2,891,245). J. H. 
Coogan and R. Klein, Arcadia, Calif., 
assignors to the U. S. Navy. 

Selective filter comprising a balanced 
modulator adapted to receive a Doppler 
signal subject to constant variations in 
frequency within a predetermined range. 
The output of a local oscillator is connected 
to the modulator, suppressing the local 
oscillator contents of its output, and pro- 
ducing a constant intermediate frequency 
signal. 

External flow jet flap (2,891,740). 
Campbell, Warwick, Va. 

Trailing edge flap on a wing rearward 

of a jet engine suspended in a pod. The 
exhaust of the engine travels upward and 
strikes the flap when the flap is in the 
deflected position. 
Self-cooled turbine drive (2,893,204). 
R. J. Anderson and J. T. Hamrick, Park- 
view, Ohio, assignors to Thompson Ramo 
Wooldridge, Inc. 


RE Exhaust Air Overboard 


Integrated cooling system for the 
generator and other temperature critical 
components. Ram air or regeneratively 
cooled turbine discharge air is circulated 
around the components within a shroud. 
Shaped cavity explosive charge (2,892,- 
407). N. A. MacLeod, Altadena, Calif. 

Cylindrical outer shell with an inwardly 
extending metal-conical liner at one end 
and a flat base at the other end. Detona- 
tion waves in a primer are transmitted to 
the main explosive charge, collapsing the 
liner progressively, and increasing the ve- 
locity of the penetrating jet. 
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Nozzle (2,894,692). J. H. Ledbetter, 
H. J. Kircher and R. J. Novotny, Rock- 
away, N. J., assignors to Thiokol Chem- 


ical 
Variable area nozzle for rockets having 


provision for controlling the gas pressures 
evolved from burning solid propellants 
in accordance with temperature. The 
nozzle throat area is varied in accordance 
with ambient temperature. 
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Elements of Physical Metallurgy, by A. 
G. Guy, Addison-Wesley Publishing Co., 
Inc., Reading, Mass., 1959, 544 pp. 
$9.50. 

Reviewed by Donatp S. CLark 
California Institute of Technology 


A review of any book is nothing more 
than a person’s opinion. There cannot 
be an absolute evaluation of a book. 
However, a reviewer can judge the logic 
of the presentation, its accuracy and its 
completeness in relation to the title. 
There are problems of evaluation even on 
this basis, since no two people may 
envision the same scope of a subject or 
the same order of presentation. This 
book presents accurate information on the 
subject of physical metallurgy. It pre- 
sents somewhat more than the elements. 
Those who master the material contained 
in this book should have a good grasp 
of the principles of the subject. 

This reviewer cannot agree with the 
order of presentation of the material. 
The quantitative approach to the subject 
is good. However, the mathematical 
approach is overdone in some instances. 
The author develops some of the relations 
in kinetic theory, but leaves the matter 
hanging in the air without direct applica- 
tion. There is an indication that the 
relation may be referred to later in a 
discussion of nucleation. It might be 
preferable to develop these relations if, 
as and when they are to be applied. 

The sequence of presentation is some- 
what confusing at times. For example, 
the chapter on phase diagrams includes 
radiography, fractures and nondestructive 
testing. Why these subjects should be 
included with phase diagrams is not clear. 
This approach tends to create a hodge- 
podge impression and may cause some 
confusion to students. 

Proper care in logical presentation is 
not maintained. The first mention of 
steel appears on page 19 but it is not 
defined at this point. The definition does 
not come until page 213, and this is 
not listed in the index. There is a tend- 
ency to mention subjects before they 
are explained. For example, reference 
to the eutectoid reaction for the hardening 
of alloys in the iron-iron carbide system 
is referred to on page 199, but the subject 
is discussed in detail beginning on page 
465. The first discussion is supposedly 
in connection with the general treatment 
of phase diagrams, whereas the latter one 
is specifically devoted to the heat treat- 
ment of steel. For some individuals this 
approach may be good; for others it may 
be a little confusing. 

The illustrations are in general very good 
and appropriate. The index has some 
deficiencies of incompleteness. The pub- 
lisher has made an error in the index in 
that “Sp’’ appears after ‘“Su.’”’ This 
causes a little difficulty in referring to 
the index. 
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Ali Bulent Cambel, Northwestern University, Associate Editor 


In general, the book is an interesting 
approach to the subject and will be found 
to be a reliable text. 


High Temperature Effects in Aircraft 
Structures, edited by Nicholas John 
Hoff, Pergamon Press, New York, 
1958, 357 + vii pp. $12. 

Reviewed by James J. KauzLaricu 
Worcester Polytechnic Institute 


This book is AGARDograph 28 pub- 
lished for and on behalf of the Advisory 
Group for Aeronautical Research and 
Development, North Atlantic Treaty 
Organization. Its purpose is to bring 
up to date the general state of knowledge 
concerning aerodynamic heating and its 
effects on aircraft structures. In accom- 
plishing this end, the editor has persuaded 
17 outstanding scientists to write review 
articles on special topics allied with the 
book’s title. (A list of topics and authors 
is given at the end of this review.) 

The several authors have not developed 
their subject in a textbook style. Rather, 
they have presented results in as succinct 
a form as possible, but being careful to 
reference all information. 

The presentation is far from wordy 
and must be read carefully. In light of 
the difficulty of reading this book, the 
publishers should have used a larger print. 

The chapter on materials by Pol Duwez 
sets the tone of the materials section of 
the book. The introduction defines three 
fundamental types of crystalline inorganic 
materials according to their atomic bonds 
as metallic, ionic and hard metals. The 
important physical properties pertinent 
to the subject are next discussed, and the 
remainder of the chapter is used to discuss 
particular materials with respect to charac- 
teristics of importance to the subject. 
The succeeding chapters on plastics and 
glasslike materials are organized in similar 
manner. 

A discussion of creep and relaxation in 
metals and of plastics follows. The 
problem of creep is considered from a 
mechanical model point of view, where it 
is shown that an arrangement of springs 
and dashpots can describe the behavior 
of actual solids for a sequence of stresses 
or strains. Much space is used in ex- 
plaining the form of the creep curves in 
terms of dislocation theory of the plastic 
deformation of crystals. The discussion 
proceeds to a review of the strengths that 
may be exhibited by available materials 
under typical conditions. 

Another important section of the book 
is concerned with thermal buckling and 
creep buckling. The various theoretical 
approaches are discussed and applied to 
idealized problems. The analysis is com- 
pared to experimental results, and it is 
brought out that much research is needed 
on this subject. 

This reviewer cannot help but wonder if 
this book is really making an important 


contribution to the scientific literature 
Actually there is a great deal of simple 
cataloguing of information. In some sec- 
tions there is a strong indication of the 
author’s insight into his subject, but th 
abbreviated form of the presentation doe 
not foster much of this. Of course 
almost every chapter could easily have 
been expanded to book size. 

The complete listing of chapters witl 
respective authors are: 1. Introduction 
N. J. Hoff, Stanford University; 2. Ex- 
ternal Sources of Heat, Martin H. Bloom, 
Polytechnic Institute of Brooklyn; 3 
Heat Transmission in the Structure 
Frederick V. Pohle, Polytechnic Institute 
of Brooklyn; 4. Materials for High Tem- 
perature Aircraft Structures, Pol Duwez, 
California Institute of Technology; 5. 
Non-metallic Structural Materials at High 
Temperatures, C. Gurney, University 
College of South Wales and Monmouth- 
shire; 6. Glass-like Structural Materials at 
High Temperatures, Ivan Peyches, Cie. St. 
Gobain; 7. Creep and Relaxation in 
Metals, N. P. Allen, National Physical 
Laboratory; 8. Creep and Stress Relaxa- 
tion of Plastics, William N. Findley, 
Brown University; 9. Fatigue of Struc- 
tural Materials at High Temperatures, 
B. J. Lazan, University of Minnesota 
10. Thermal Stress, L. Broglio and P. 
Santini, University of Rome; 11. Buckling 
Caused by Thermal Stresses, A. van der 
Neut, Technische Hogeschool Delft; 12. 
Stress Distribution in the Presence of 
Creep, N. J. Hoff, Stanford University; 
13. Creep Buckling, B. Fraejis de Veubeke, 
Université de Liége; 14. Influence of 
Aerodynamic Heating on Aeroelastic 
Phenomena, R. L. Bisplinghoff and John 
Dugundji, Massachusetts Institute of 
Technology; 15. Experimental Methods 
in High Temperature Structural Research, 
J. Taylor, Royal Aircraft Establishment; 
16. Models and Analogs, R. R. Helden- 
fels, NACA Langley Aeronautical Labora- 
tory. 


Guided Missiles, Operations, Design, and 
Theory, sponsored by the Dept. of the 
Air Force, McGraw-Hill Book Co., 
Inc., New York, 1958, 575 pp. 

Reviewed by LAWRENCE S. Brown 
Grumman Aircraft Engineering Corp. 


This book is a reproduction of an Air 
Force manual on guided=missiles funda- 
mentals, compiled in( 1955 by the instruc- 
tors and staff of the Air Training Com- 
mand’s school for missiles at Lowry 
A.F.B., Colorado. In the foreword, Lt. 
Gen. C. T. Myers, USAF, states: ‘While 
seemingly complex, this book is little 
more than a primer in a field where science 
and technology move in their highest 
form.’ This appears to be a fair state- 
ment although ‘‘voluminous”’ would be a 
better word than ‘“complex.’’ 

Coverage of the subject is quite com- 
plete as indicated by the following para- 
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phrase of the contents: The story of 
guided missiles, aerodynamics, propulsion, 
physics, control components, guidance 
components, control systems, guidance 
systems, trajectory considerations, tactics 
and instrumentation. The book is of 
high potential value to the general tech- 
nician or to the specialist engineer who 
wants a briefing (or a reference) on matters 
other than his specialty. It is unfortunate 
that the first chapter was not de-oriented 
from Air Force usage. In the recount 
of foreign developments, the name of 
Wernher von Braun is conspicuously 
omitted. The last half of the chapter 
describes the USAF organization for de- 
velopment of guided missiles, but is silent 
about how the administration of develop- 
ment of a Polaris or Pershing is handled. 
This is, of course, of little scholastic im- 
portance, but tends to annoy an unbiased 
taxpayer. 

The technical treatment of the balance 
of the subjects is excellent, although not 
uniform in level of coverage. The section 
entitled: ‘How optical and electronic 
principles work together’ (following a 
section on Optics) is particularly unique, 
as is the section on computer units of 
guidance systems. The illustrations— 
schematics, block and circuit diagrams, 
and pictorial—are far superior to those 
of most similar publications. 


Progress in Elementary Particle and 
Cosmic Ray Physics, Vol. IV, edited by 
J. G. Wilson and S. A. Wouthuysen, 
Interscience Publishers Inc., New York, 
1958, 470 pp. $12.50. 

Reviewed by S. F. Sincer 
University of Maryland 


This is the fourth volume in a series 
which deals with cosmic ray physics. 
However, the subject is now divided into 
two main branches. One is concerned 
with the nuclear interactions of the ex- 
tremely high particles which the cosmic 
radiation furnishes us, and with the new 
and strange elementary (i.e., subatomic) 
particles produced in these nuclear inter- 
actions. Thanks to the development of 
high energy accelerators some of these 
nuclear interactions can now be repro- 
duced in the laboratory, and it was 
therefore thought best to discuss the 
elementary particle physics using both 
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Propulsion and Power 
(Non-Combustion) 


An Ion Source Operated by High Fre- 
quency Pulses, by G. S. Malkiel’ and 


Eprror’s Note: Contributions from Pro- 
fessors E. R. G. Eckert, J. P. Hartnett, T. 
F. Irvine Jr. and P. J. Schneider of the 
Heat Transfer Laboratory, University of 
Minnesota, are gratefully acknowledged. 


136 


cosmic ray and accelerator data. The 
first three chapters are concerned with 
this general area. The first by D’Espagnat 
and Prentki of the European Nuclear 
Center in Geneva discusses’ aspects of the 
strong interactions of the new particles. 
This is a highly technical paper which 
reviews the existing theories. Strong 
interactions are of the type that occur 
in connection with nuclear forces, whereas 
the so-called weak interactions are as- 
sociated with radioactive data decay. 

The second chapter by Walker of the 
University of Wisconsin discusses prop- 
erties and production of K-mesons. The 
K-meson is a particle discovered only 
about 10 years ago; its mass is approxi- 
mately 900 times that of the electron. 
It decays into lighter mesons, the so-called 
pi-mesons and mu-mesons which in turn 
decay into electrons and neutrinos. Be- 
cause of the short lifetime of the K- 
particles (of the order of 10~* sec) their 
study is extremely difficult. There exist 
also neutral K-mesons which have a life- 
time of 10~" sec. 

The third chapter by Fowler and Wolf- 
endale of the University of Bristol and 
the University of Durham, England, re- 
spectively, deals with interactions of 
minu-mesons which were the first mesons 
discovered. Their interaction, however, 
with nuclei is very weak, and for this 
reason the number of interactions is small 
and experiments are rather difficult. 
However, when a mu-meson is captured 
by a positive nucleus into an orbit, then 
the chance of a nuclear interaction be- 
comes greater because of the close prox- 
imity and the long time available. In 
many respects, therefore, a mu-meson 
behaves like a heavy electron, its mass 
roughly 210 times that of the electron. 

The next two chapters deal with an 
aspect of cosmic rays which is becoming 
of importance to astrophysicists. In 
studying the time variations and origin of 
cosmic rays we really deal with the problem 
of how these highly energetic particles are 
produced in nature, in other words, the 
electromagnetic properties of the solar sys- 
tem and the galaxy. Areview of geomag- 
netics is first, since it forms the basis of 
much of the analysis of the cosmic radia- 
tion. The first half of the review deals with 
our present knowledge of the primary 
radiation, its composition the energy 


spectra of the primary protons, the 
primary alpha particles and the heavier 
primaries. An outstanding problem with 
the primary radiation is the occasional 
absence of low energy cosmic rays, and 
one of the puzzling features, of course. 
is the mechanism which removes them 
from the solar system, or at least from 
the vicinity of the Earth. The second 
part discusses the time variations of 
cosmic rays, in particular the increases 
associated with solar flares, the decreases 
associated with some magnetic storms and 
variations during the sunspot cycle. 
Attempts are made to account for the 
existing theories and to present a co- 
herent picture of the reason for these 
time variations. It turns out that these 
time variations give us a tool of studying 
the interplanetary space, in particular 
the existence of magnetic and electric 
fields in the interplanetary space. The 
cosmic rays, being charged particles, are 
deflected by magnetic fields and thus act 
as probes. Of particular interest to space 
operations are the discussions of the in- 
tensity increase with altitude, its depend- 
ence on latitude, the penetration charac- 
teristics of primary cosmic rays as a 
function of their energy and charge, and a 
discussion of the cosmic ray albedo. 

The last chapter is by Ginzburg of the 
USSR Academy of Science in Moscow. 
It deals with the origin of cosmic radiation 
and describes mainly the theoretical work 
of Shklovski and of Ginzburg in trying to 
tie together the production of cosmic rays 
and production of radio noise in one 
theory. These views have been strikingly 
confirmed in recent astronomical work on 
the polarization of the light from the Crab 
nebula supernova. The chapter discusses 
the motion of cosmic ray particles in inter- 
stellar space and their interactions with 
interstellar gas. The electrons which 
are the final end product of these oc- 
casional nuclear interactions are then 
held to be responsible for the emission of 
radio noise by gyrating in the interstellar 
magnetic fields. Finally, the possibility 


that supernovae and novae are the actual 
sources of cosmic rays are discussed. It 
is interesting to see that many of the out- 
standing problems defined by the authors 
have found their experimental realization 
in the cosmic ray 
Sputnik ITT. 
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1958, p. 751. 

On Formation of Magneto-Hydrody- 
namic Shock Waves, by 8. Segre, Nuovo 
Cimento, vol. 9, no. 6, Sept. 16, 1958, pp. 
1054-1057. 


Viscous Dissipation in Low Prandtl 
Number Boundary-Layer Flow, by E. M. 
Sparrow and J. L. Gregg, J. Aero/Space 
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the Study of Leading Edge Flows, by J. A. 
Laurmann, Phys. Fluids, vol. 1, no. 6, 
Nov.-Dec. 1958, p. 469. 
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367. 

Supersonic Profiles with Minimum 
Drag, by Iu. D. Shmyglevskii, J. Appl. 
Math. Mech. (Translation of Prikladnaya 
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(Translation of Prikladnaya Matematika i 
Mekhanika), vol. 22, no. 2, 1958, 375- 
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On the Motion of Thin Rigid Bodies 
Under the Influence of Shock Wave 
Forces, by S. S. Grigorian, Akad. Nauk, 
SSSR, Izvestiia, Otdelenie Tekh. Nauk, 
Mekhanika i Mashinostroenie, no. 1, 1959, 
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On Ablation for the Recovery of Satel- 

_ lites, by S. Georgiev, H. Hidalgo and 

-M. Adams, Avco Mfg. Corp., Avco Res. 
- ‘Lab., Res. Rep. 47, March 1959, 14 pp. 


Flight Mechanics 


On the Solution of a Degenerate Varia- 
tional Problem and the Optimum Climb of 
a Cosmic (Space) Rocket, by V. A. Egerov, 
J. Appl. Math. Mech. (Translation of 
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